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Both total and diffractive cross sections from HERA are successfully confronted with JIMWLK 
evolution equations in the asymptotic pseudo-scaling region. We present a consistent, simul- 
{S) taneous description of both types of cross sections that includes NLO corrections in the form 

of running coupling and energy conservation corrections. The inclusion of energy conservation 
^ corrections allows to match all available data with x^j < -02 i.e. up to Q 2 < 1200 GeV 2 . We 

discuss the effects of quark masses including charm, contrast asymptotic and pre-asymptotic fit 
strategies, and survey non-perturbative uncertainties related to impact parameter dependence. 
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£ 1 Introduction 

Much of the abundant particle production in modern collider experiments at high energies is frig- 
id gered by gluon channels, thus imprinting the features of gluon phase space on many observables. 

This is the basis of the importance of the Color Glass Condensate (CGC) [1-29] for virtually all 
current collider experiments be they designed to answer particle- or and heavy-ion-physics questions. 

The characteristic feature of enhanced gluon emission into the final state at high energies is the 
emergence of an energy dependent transverse correlation length R s (x). Its associated conjugate 
momentum scale Q s (x) ~ 1/R s (x) signifies the onset of gluon saturation, hence the name saturation 
t-H scale. As gluon numbers rise with energy, the correlation length of the dense gluon cloud shrinks, 

the saturation scale increases. 
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The mere existence of such an energy driven scale has led to a large body of phenomcnological 
literature, often based on identifying quantities that can be expected to crucially depend on the 
saturation scale Q s to gain insight into the energy dependence of some observable by applying a 
scaling argument. The origin of this idea predates the observation that perturbative QCD allows 
us to predict Q s -scaling in the context of the JIMWLK equation (or in an independent scattering 
approximation the BK equation [21-25]) and has provided us with the discovery of what was called 
geometric scaling in HERA data [30-32]: plotting the ep cross sections measured at HERA not 
as a function of rapidity Y = ln(l/x) and momentum transfer Q 2 independently, but instead 
as a function of the scaling variable Q 2 /Q 2 (x) reveals beautiful scaling features of the data for 
x < 10~ 2 that extend even to diffractive measurements. This scaling subsumes the complete 
energy dependence of the data at x < 10 -2 in a single energy dependent function, Q s (x). 
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The standard approach to cross correlate HERA data in the Q 2 direction using Dokshitzer-Gribov- 
Lipatov-Altarelli-Parisi (DGLAP) evolution equations reproduces this scaling feature only acciden- 
tally, in the sense that they match data which themselves exhibit scaling. There is no intrinsic 
reason why such a scaling feature should emerge within the domain of validity of the linear DGLAP 
formalism 1 . By contrast, in the CGC context, scaling is a natural byproduct of nonlinear features 
of gluon emission and saturation, which adds a particular interest to a comparison of CGC results 
with HERA data. 

The derivation of the JIMWLK evolution equation relies on the prerequisite that there exists a 
frame in which the gluon field of the target becomes strong, while the gluon field of the projectile is 
weak. The resummation of the large target field then induces nonlinearities that capture the effects 
of gluon saturation and, through evolution, driven by perturbative gluon emission in the weak field 
projectile, induce saturation and the appearance of a saturation scale Q s (x) (its inverse R s ~ 1/Q S 
has the interpretation of a transverse correlation length in the dense gluon cloud) . 

While the mere presence of such an energy dependent scale will impose its mark on any observable, 
strict mathematical scaling, i.e. the notion that all the cross sections and correlators of the theory 
share an energy dependence that can solely be expressed in term of that of the saturation scale 
Q s (x) is too naive an expectation: It is more natural that different n-point functions have their own 
associated scales unless tied together by JIMWLK evolution by belonging to the same Balitsky- 
hierarchy for which a strong factorization feature holds that allows to express these generically 
independent correlators to be expressed in terms of a single one, at least to good approximation. 
The scaling observed in the HERA total cross section should be interpreted as a signal for such a 
correlator factorization and the possibility to indeed truncate the Balitsky hierarchy of the dipole 
amplitude, the correlator dominating the total cross section. 

Since the same correlators dominate both total and diffractive cross sections at HERA, the same x- 
dependent scale should also be in effect there, a fact established early in the Golec-Biernat-Wiisthoff 
(GB-W) model [30-32]. 

A second caveat arises from the QCD scale anomaly: even for observables which are predominantly 
determined by only a single correlator, true scaling is only to be expected at leading order (LO). 
At next to leading order (NLO) exact scaling receives small corrections [which manifest themselves 
as a slow drift of correlator shapes] purely induced by the scale anomaly, i.e. the running of the 
coupling. We will refer to this regime as the asymptotic or (pseudo-) scaling regime. 

Despite the phenomenological success of models that assume (pseudo-) scaling of the cross section, 
it is by no means clear that the data would show clear evidence of strict or near scaling in HERA 
data or even for the HERA total cross section alone. 

The question if the scaling observed in the 7*p cross section [34-36] and the rapidity gap events [37 
40] (see [41] for the most recent combined update of all ZEUS data) at HERA affects all corre- 
lators, i.e. represents true (pseudo-) scaling, or only signifies the presence of an intrinsic scale, 
indicative merely of the presence of nonlinearities, with scaling on a correlator level only apparent 
and potentially limited to a very narrow kinematic range is surprisingly hard to settle. 

Let us emphasize that a survey of the literature presents us with a very ambiguous picture: All 
models from GB-W to the BK inspired parametrization of the dipole cross section by Iancu, Itakura 

1 The derivation of scaling "within" DGLAP in [33] skirts the region of applicability of the argument: the scaling 
solutions shown there push into the region of large gluon densities where nonlinear/higher twist corrections are bound 
to become important. 
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and Murder [42] (IIM) to the BFKL+saturation boundary model of Mueller and Triantafyllopou- 
los [43,44] (MT) exhibit either strict or pseudo-scaling of correlators, since scaling in this sense is 
one of the main constraining features in the construction of such models. 

However, the shapes of these various (pseudo-) scaling solutions are notably different from each 
other and, as we will see below (see Fig. 7), from the scaling solution imposed by the evolution 
equation at NLO: the imposed scaling shapes in these models do not resemble the scaling shapes 
that emerge as solutions of the evolution equation. 

Moreover, even within a CGC framework, strict or near scaling of correlators does not seem to be 
required to obtain a successful fit to the HERA total cross section: This was shown by Albaccte, 
Armesto, Milhano, and Salgado in [45,46]. They omit all NLO contributions beyond running 
coupling, but start evolution with a shape close to the GB-W parametrization. In this treatment, 
correlators only approach the scaling shape imposed by the evolution equation at the far end of the 
x range covered by the HERA experiments. 

Our own simulations [47-49] add in energy conservation corrections to cover NLO corrections beyond 
the running coupling contributions and provide an equally convincing fit in the pseudo-scaling 
region. These same NLO corrections restrict our treatment to the region near (pseudo-) scaling due 
to stability and self consistency considerations that only emerge once all these NLO corrections are 
taken into account (see Sec. 2.3.3). 

We will attempt to summarize the status of the theoretical tools presently available (state of the 
art are NLO evolution combined with LO impact factors, for a more in depth discussion see below) 
and explore what kind of tension available data pose on our theoretical analysis by considering both 
the HERA total cross section and the rapidity gap events. 

As our first set of results emerges from fits to the total cross section: we will argue that simulations 
that include NLO corrections in the form of running coupling corrections and energy conservation 
corrections (the most complete set of NLO corrections currently available) favor fits in the pseudo- 
scaling region based on fits of the total cross section. The treatment of quark masses has proven 
to be somewhat problematic in earlier fits, which have generically used constituent like values of 
around 140 MeV for light quarks and 1.4 GeV for charm. Conceptually, quark masses are subleading 
in the small x limit where factors of the form (a s ) n+m (hx(l/ x)) n are used to sort contributions by 
importance and hence should not pose a serious difficulty. Choosing current quark masses instead 
of constituent quark masses we find that fits are indeed feasible with the quark masses having 
their largest effect in the nonperturbative range with Q 2 < 1 GeV 2 . We find that the idea of [30] 
to address the situation by replacing x by some x c s to modify the small Q 2 limit can improve 
the fit, although the specific form introduced in [30] proves unusable. Such resummations are by 
construction nonperturbative and should at this stage be taken to merely indicate the relevance of 
nonperturbative input in this range of phase space. 

A second set of insights emerge from a study of rapidity gap events: Our analysis is based on 
the fit parameters extracted from total cross sections and clearly shows us the limits of a fit with 
incomplete NLO input, despite the quite satisfactory fit quality. A comparison with data clearly 
requires a ggg-component in the impact factor for which at present we only have a rough substitute 
which is only reliable in the large Q 2 limit and was already devised long ago in the context of 
the GB-W model [31]. This prevents us from using observables beyond the total cross section to 
get a closer look at the details of JIMWLK evolution - this crude treatment is not suitable for 
such a precision study. It turns out that non-perturbative aspects which enter through the impact 
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parameter dependence adds additional uncertainties as one steps beyond the total cross section. 
It affects the relative normalization of individual Fock-space components as soon as NLO impact 
factors start to play an important role in resolving the structure of the cross section. 

This leads us to conclude that to address questions such as precision fits on initial conditions, 
the reliability of truncations of the full JIMWLK framework, or the size and nature of sublead- 
ing iV c -corrections as advocated in [50] with any definiteness, we need full knowledge of all NLO 
contributions including the appropriate impact factors for the observable in question as well as an 
improved understanding of the non-perturbative aspects of the impact parameter dependence. 

Any progress in this respect will also improve the utility of HERA fits as an input for fits to RHIC 
and LHC data and will be one of the most important tasks for the near future. 

The structure of the paper is as follows: We begin by recapitulating the theoretical ingredients 
necessary to define the underlying observables at zeroth order in Sec. 2.1 and 2.2, putting some 
emphasis on the approximations underlying the expressions usually given in the literature. We 
review our present knowledge of the NLO corrections and how to use truncations of JIMWLK 
evolution to efficiently implement them in Sec. 2.3. 

Sec. 3 is devoted to the systematics of a fit to the total cross section in the asymptotic pseudo- 
scaling region. We discuss general features of the asymptotic fit such as evolution speeds, Q s in 
HERA phase space and correlator properties in Sec. 3.1. This is followed by a thorough study of 
the role of the energy conservation correction (Sec. 3.2), and the effect of quark masses (Sec. 3.3). 

The fit obtained in Sec. 3 is then applied to diffractive data in Sec. 4. We begin our discussion in 
Sec. 4.2 with estimates of the non-perturbative uncertainties induced by our lack of knowledge of 
impact parameter- (b-) dependence of the eikonal correlators in a proton or nuclear target. 

In Sec. 5 we compare our asymptotic approach that includes both running coupling corrections and 
the energy conservation correction with more conventional fits that only use the running coupling 
effect, but leave out the energy conservation correction and use the pre- asymptotic regime of small x 
evolution. We show that features of the solutions that are recovered perturbatively in the asymptotic 
region must be imprinted at least partially via the initial condition in the pre-asymptotic approach. 

Sec. 5 collects our main results in attempt to provide a synthesis. 

Several appendices provide a number of generic expressions (Sec. A), ancillary results partly in- 
dispensable to reconstruct our numerical simulations (Sec. B through D), as well as number of 
consistency checks (Sec. E). 

2 Deep inelastic scattering at small x 
2.1 Total 7M cross sections at zeroth order 

There are many phcnomcnological applications to various physical processes and differential cross 
sections that are based on the idea of obtaining energy dependence from the scaling behavior in 
terms of the saturation scale. However, the observable directly addressed by JIMWLK evolution is 
simply the total cross section in very asymmetric collisions such as eA or pA experiments. 

The strong asymmetry between projectile and target serves to justify the notion that the gluon 
field of the nuclear target with atomic number A can be thought of as much larger than that of 
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the much simpler electron or proton projectile. The asymmetry is used to describe the projectile in 
terms of a very simple wave function with only a very few (valence) partons, which then scatter on 
the large target field with a very large longitudinal momentum. This justifies the description of the 
interaction of such a projectile constituent with the large target field in a no recoil approximation. 
The no recoil approximation fixes the projectile constituent onto a worldline at a transverse position 
unaltered during the interaction with the target. Multiple interactions with the target field then 
build a non-Abelian eikonal factor, a path ordered exponential that captures the interaction of the 
constituent with the target field. This way the interaction of a quark in the projectile with the target 
is represented as an SU (N c )-vahied field in the transverse plane U x , an antiquark analogously enters 
as a U%. and a gluon in the projectile interacts as U x (the tilde denotes the adjoint representation). 
The graphical notation used to represent this is shown inf Fig. 1 for the example of the zeroth order 
contribution to the "f*A scattering amplitude. 
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1: Diagrammatic representation of the amplitude for j*A scattering at small x at momentum transfer 
— q 2 . Light cone "time" x~ runs from right to left. The interacting "out-state" (left) contain nontrivial 
interaction between projectile and target, which is marked by a vertical bar (blue online) at i~ =0 that 
indicates the interaction region and markers for the Wilson lines picked up by the projectile constituents. The 
non-interacting "in-state" (right) instead has no interactions and correspondingly trivial Wilson line factors at 
x~ = 0. 



The corresponding zeroth order total cross section arises as the absolute value squared of the 
difference of this interacting diagram with its noninteracting counterpart 2 
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(2.1a) 



If the transverse momentum integrals are unrestricted, they will identify the transverse coordinates 
left and right of the cut, so that the the [/-content of these diagrams is partially simplified: 



(2.1a) = ' j j 

ttU y , U^, Ul m ~ >tc > trl 



+ 



(2.1b) 



trUvUX 



2 Note that we keep the dashed vertical line that marks x =0 also in the noninteracting case. We will need this 
below to distinguish where gluon vertices connect with respect to x~ = 0. 
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Note in particular that, with this assumption, the {/-factors left and right of the cut in the out-out 
overlap cancel against each other. Diagrammatically, we have 

(2.2) 

While this is assumed universally in the literature without discussion, one should be clear that 
this is an approximation: phase space integrals over the transverse momenta ki of the final state 
quarks are limited by W 2 :— s 7 *a (the invariant mass of the 7*A-system) both in t := — (fei + k^) 2 
and the invariant mass of the produced pair. Consequently, also the momentum integrals over the 
transverse momenta are not unrestricted and thus the primed and unprimed coordinates are not 
quite the same in an exact treatment: even the total cross section will have a contribution from a 
four Wilson line operator. At small Q 2 /W 2 , and presently achievable theoretical accuracy (see our 
discussion of NLO contributions below and the discussion of coincidence limits given in [50,51]), 
it is fully justified to follow custom and take the limit shown in (2.2) and to incorporate this four 
point function only in a coincidence limit in which it becomes trivial. 

With this caveat, it is the {/-content of the two remaining diagrams, namely 
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S% := \ T VJ (2.3) 



and its complex conjugate that capture the interaction of the qq pair with the target. 

The presence of the target wave function induces an energy dependent averaging process, that for 
the applications below will be strictly real, so that we may (and will) not distinguish Sy(x, y) := 
(Sxy)(X) from its complex conjugate {S y q x ){Y) in the following. For the total cross section, the 
target interaction of Eq. (2.1) can be fully summarized by the dipole operator 

Nil ■= w c tl{1 ~ UxU y ] (2 ' 4) 

(and its complex conjugate), and the average Ny{x,y) :— (N^ y ) (Y) , the dipole amplitude. It 
is the average over the target wave function, an operation that involves both perturbative and 
non-pcrturbativc information, that proves the most difficult part of this calculation and induces 
the energy (or Y-) dependence of the cross section. The tool to extract this energy dependence 
is the JIMWLK equation, and its associated framework. At zeroth and leading order (LO) in 
a s ln(l/x) this allows us to describe the 7* A cross section at a given energy entirely in terms of this 
simple amplitude without reference of "higher" Fock-space components of the projectile, simply 
by subsuming (in the sense of a renormalization group procedure) all other strongly interacting 
components into the averaging procedure. 

The last ingredient of Eq. (2.1) not yet spelled out analytically, the wave function of the virtual 
photon, is known exactly [see for example [ ] and Eqs. (2.6)], with both longitudinal and trans- 
verse polarizations contributing additively to the total cross section, a^ ^(Y,Q 2 ) = ctJ p (Y,Q 2 ) + 
o-T P (Y,Q 2 ). 

This leads to an expression for the cross section in the form of a convolution in terms of the 
transverse coordinates which characterize the eikonal scattering position of the qq pair. Using 
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r = x — y and b = zx + zy to denote dipole size and impact parameter respectively, one obtains 
(notations inspired by [53]) 

1 

ai P L (Y,Q 2 ) = J2 \ Sr Jdz^ L (z,r,r\Q 2 ) J d 2 b2Nf(x,y) , (2.5) 
f o 

where z and z := 1 — z denote the longitudinal momentum fractions carried by the quark and 
antiquark respectively. For a fixed Q 2 , polarization and flavor /, the photon wave function product 
L (z, r, r' , Q 2 ) 3 encodes the probability to find a qq pair of size |r|, polarization T or L, and 
longitudinal momentum fraction z inside the virtual photon: 



& T ( Z) r,r';Q 2 ) = ^yej^ + z*)Qj Kl {Q ^^{d } \r' \) 

+ m)K {Q f \r\)K Q {Q f \r l \)) , (2.6a) 

${(z,r,r';Q 2 ) = ^l e 2 4Q 2 z 2 z 2 K (Q f \r\)K (Q f \r'\) . (2.6b) 

In the above, et and m/ denote the charge and mass of the quark with flavor / and 

Q 2 f = zzQ 2 +m 2 f . (2.7) 

The impact parameter (b) integrated dipole amplitude has the interpretation of a qq-dipole cross 
section on the target: 

a q<] (Y, (x - y) 2 ) := 2 J d 2 b Nf{x, y) . (2.8) 

It carries the energy dependence of the cross section in terms of 1/x — e Y , the relative boost 
factor between the projectile and the target. The separation into wave-function factors (generically 
called impact factors) and dipole cross section (more generically Wilson line n-point functions) is 
prototypical to all observablcs in the high energy limit and extends to higher order in perturbation 
theory. 

A description in terms of structure functions, F%, Ft^l corresponds to a purely kinematical reparametriza- 
tion according to the standard relation 

aTj{x hh Q 2 ) = A -^^F 2 {x hh Q 2 ) = (F T (* bj ,Q 2 ) +F L ( Xhj ,Q 2 )) . (2.9) 

Please note that this does not imply a general link with particle distributions outside the parton 
gas region with Q 2 3> Q 2 {x) where a twist expansion becomes valid. 



3 The notation for the expressions shown here and in Eq. (2.6) anticipate the diffractive case. There the structure 
of the wave function overlaps remains the same, only the dipole sizes on both sides of the final state cut must be 
distinguished. 
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2.2 Rapidity gaps in 7* A at zeroth order 



Apart from the total cross section, large complementary data set is available for rapidity gap 
events, in which the virtual photon fragments into (predominantly) a qq-pair accompanied by a 
gluon shower (which then hadronizes before it reaches the detector) that remains well separated 
from the target fragmentation region by a large rapidity gap. The kinematical setting for rapidity 
gap events is sketched in Fig. 2. 
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Fig. 2: Rapidity gap events differ from generic events contributing to the total cross section by a target side 
rapidity gap of size Y gap = ln(l/a;p) into which no gluons are emitted. This gap is complemented by a projectile 
fragmentation range of size Yf rag = ln(l//3), such that Y = Y gap + lf rag . 



The experimental situation is characterized by a quite strong similarity of the energy dependence 
of both total and diffractive cross sections. 

Theoretically, all differences with the expression for the total cross section Eq. (2.1) arise from 
restrictions on the final state: With a rapidity gap on the target side, the target stays intact. 
Therefore, in the final state, the target is projected back onto its wave-function in each amplitude 
factor. This leads to separate target averages (...)(Y) in each amplitude. Since no net color 
can be exchanged across the gap, the qq-Rn&l state is necessarily projected onto a singlet. Once 
perturbative corrections are taken into account, and partons added to the projectile fragmentation 
region, they remain in an overall singlet. 

At zeroth order in a s ln(l/;r), the projectile only contains a qq pair, the cross section is given by 

(^T--T)^- C (T--T) 

in close analogy with (2.1). 

We note that each of these restrictions individually (the separate averages as well as the singlet 
projection of the projectile constituents) will prevent the simplification of the [/-content that takes 
place in the out-out overlap for the total cross section due to Eq. (2.2). Instead 



As an immediate consequence, not only the mixed (in-out) overlaps but also the out-out overlap 
will acquire nontrivial energy dependence. 

What is left to understand are the restrictions on the transverse coordinates of Wilson lines in this 
expression as imposed by phase space integrals - again such a restriction alone would be sufficient 
to induce a nontrivial Wilson-line four-point function. To this end, note that j3 and Q 2 together 
determine the invariant mass of the projectile fragments: 

O 2 

P=^r^ o- (2-12) 

This restricts the integration over transverse momenta. To be specific, assume n projectile con- 

n 

stituents, denote longitudinal momentum fractions by Zi (with ^ Zi = 1) and final state transverse 

i 

n 

momenta by fc; := ki + z^A (where A :~ ^ki). Then the invariant mass of the n-particle final 

i 

state is 

« £.2 , 2 

Z>i 

I 

so that a restriction on Mx imposes one linear constraint on the ft?. To understand what this 
implies for the transverse coordinates, define 

T<i . - Xi — X n ^ .— X^ 3? n , b . ^ ^ ZjT*2, b . — ^ ^ Zj^T , (2.14) 

i i 

and rewrite the exponents of the transverse momentum phase factors as 

n n—1 

E k i ■ ( x i - *i) = A • (b' - b) + E ^ • (r'i - n) • (2.15) 

i i—1 

This implies that integration over momentum transfer t — —A 2 will identify the light cone cm. 
coordinates b and b' if one ignores the kinematical upper limit on t as for the total cross section. 
One of the n — 1 ki integrations on the other hand is restricted by (2.13), so that one of the n—1 
independent distance pairings t\ and r[ will remain independent after all unconstrained integrals 
are carried out. For the qq final state encountered at leading order, there is only one independent 
momentum variable available to begin with k = k\ = — ki- It is, therefore, tied directly to Mx via 

Ml, = (2.16) 

and leaves behind phase factor containing (r' — r) ■ k where the dipole sizes r and r' in amplitude 
and complex conjugate amplitude remain independent. For each flavor independently, the length 
of k is fixed in terms of (3 and the quark mass via (2.16) and (2.12) as 

K 2 f := zzQ 2 ^- - m 2 f . (2.17) 
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One obtains (see again [53]) 



d3 




(2.18) 



where Y = ln(l/x). Integrating this result over /3 leads to the total diffractive cross section. This 
step identifies r with r' if one extends the upper phase space boundary in Mx from W 2 to oo as 
discussed earlier. The result in turn maps back onto the total cross section if one removes the singlet 
projection in the final state which replaces (Nf(r', b'))*N*?(r, b) by (N™(r',b'))* + Nf(r, b). 

Note that this relationship of cross sections uniquely identifies the x arguments of the averages in 
the diffractive cross sections to be the overall Bjorken x of the process, not x-p as usually assumed 
in the literature. The numerical effect of such a replacement is, however, not large enough to affect 
the quality of any diffractive fits with state of the art expressions. These expressions suffer from 
more serious defects: incomplete NLO impact factors and nonperturbative normalization effects 
associated with the corresponding impact parameter averages as will be discusses below. 

2.3 Beyond zeroth order 

Leading order (LO) corrections to the above resum contributions proportional to (a s ln(l/x)) n and 
are fully taken into account by solving the LO-JIMWLK equation. The impact factors receive no 
corrections to their zeroth order form. At next to leading order, when contributions proportional 
to a s (a s \n(l/x)) n are taken into account, both JIMWLK evolution and the impact factors receive 
corrections. 

Running coupling corrections to the evolution of Wilson line n-point functions have been calculated 
in full generality [54-56]. The remaining conformal corrections to evolution have been obtained 
by Balitsky and Chirilli [ ], who presently work on the expression for NLO impact factors. A 
generalization of both aspects for arbitrary Wilson line n-point functions is yet to be devised. These 
ingredients would be required to extend the only existing treatment of JIMWLK-evolution [58] 
beyond leading order. 

Fortunately truncations of JIMWLK evolution to finite sets of evolution equations, such as the 
Balitsky-Kovchegov (BK) equation or its more general GT (Gaussian truncation) counterpart allow 
us to use the available information to implement evolution at NLO. Accuracy on the impact factor 
side at present remains at LO. NLO accuracy for the impact factors introduces terms including 
^(/-correlators into the expressions for both the total cross section Eq. (2.9) and the diffractive 
cross section (2.18). We will see below that fits to the total cross section are not affected strongly, 
but that already the description of rapidity gap events suffers noticeably from this limitation. 

Even with more modest goals in mind, such as the description of the total cross section it is 
mandatory to include NLO effects at least on the level of evolution equations. Here NLO corrections 
induce qualitatively new effects like scale breaking and a quantitatively important reduction in 
evolution speed compared to the LO situation. With the impact factors remaining at LO it is 
sufficient to know the Independence of the dipole cross sections entering Eqs. (2.9) and (2.18), 
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but once the NLO impact factors are known this is no longer sufficient (qqg operators will require 
consistent treatment) and one is forced to either use the full JIMWLK-evolution framework or 
choose a "suitable" truncation. To appreciate what is involved, we briefly recapitulate the necessary 
tools. 



2.3.1 JIMWLK and its truncations at LO 

The JIMWLK evolution equation provides a means to calculate the energy- (or rapidity-) depen- 
dence of arbitrary U -correlators by first introducing an energy- (or rapidity-) dependent statistical 
weight Zy[U] for the configurations of the [/-fields. The dipole correlator of Eq. (2.4) is then 
expressed as a functional integral of the form 

(tr(l - U x Ul)) (Y) = j D[U]^-tv(l - U x Ul)Z Y [U] (2.19) 

where D[U] is a functional Haar measure. This is meaningful in the sense that it allows to calculate 
the average of any operators, if it is possible to describe the evolution of all averages in terms of 
the evolution of the weight Zy[U] defining the averaging procedure. This is the main content of 
the JIMWLK equation: it abstracts the energy dependence of the average from the operator being 
averaged by describing it as a functional evolution equation for Zy[U]. At LO, the equation takes 
the form of a functional Fokker-Planck equation 

JyZ Y [U] = -H, 1M wlkZ y [U] (2.20) 

that traces how additional gluons are added to the phase space of the projectile as one increases the 
energy of the collision. To arrive at (2.20) one needs to prove [18] that this equation indeed allows 
to find the energy dependence of arbitrary [/-correlators, not just the simple qq-operator entering 
the dipole cross section. For such generic operators 0[C/], Zy[U] defines a target average via 

(0[U])(Y) := J D[U]6[U}Z Y [U] (2.21) 
so that (2.20) implies an evolution equation for each such operator 0[U) that takes the form 

4^{0[U]){Y) = -(H nMWLK d[U})(Y) . (2.22) 
dY 

One of the features of JIMWLK evolution is that non-singlet correlators are exponentially sup- 
pressed by infrared divergent contributions, only singlet correlators survive. 

A key feature of ^-evolution is that Eq. (2.20) gives rise to coupled hierarchies of evolution equations 
for [/-correlators, known as Balitsky hierarchies. This already becomes manifest in the evolution 
equation for the qq operator at LO: It can be written as 



d 

dY { 



tv(U x Ul))(Y) = ^ [ d 2 z K, xzy (([U z ] ab tv(t a U x t b Ul))(Y) - C f (tv(U x Ul))(Y)) 



(2.23) 
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or, using (2.3) and the Fierz identity 



[U z ] ab 2tv(t a U x t b Ul) = ix(U w Ul) tr(C/ z C/t) - —ix(U x Ul) ( 2 .24) 



MS 



2?0"vK Y ) = ^Jd 2 z IC xzy (S xz S zy - S xy )(Y) . (2.25) 
The integral kernel in both (2.23) and (2.25) is given by [3,21] 



(x — 2/)^ 

ICxzy ■= t r^—, ^ • (2.26) 

(x - z) 2 (z - y) 2 



Eqs. (2.23) and (2.25) do not represent closed equations since the evolution of (tY{U x U y )){Y) de- 
pends on an operator with an additional gluon operator insertion made manifest by the U appearing 
in the first term on the right hand side of Eq. (2.23). The evolution equation of that new opera- 
tor, ([?7 z ] ab tr(t a U x t b U y ))(Y), in turn will involve yet one more insertion of a gluon operator [7, 
iteratively creating an infinite coupled hierarchy of evolution equations, the Balitsky hierarchy of 
the quark dipole operator (2.3) [23,24]. JIMWLK evolution summarizes the totality of all such 
hierarchies, based on any (gauge invariant) combination of multipole operators but can only be 
solved numerically [51,58,59] at considerable computational cost. 

One may note that within such hierarchies one finds numerous cross- and self-referencing patterns 
that can be exposed by looking at coincidence limits of coordinates in the operators involved. Tak- 
ing ([f/ z ] a tv(t a U x t b U y )) (Y) as an example, one finds its evolution equation linked with that of qq- 
and gg-dipole operator averages to which it reduces in the limits x — > y and z —> x (or y) respec- 
tively. More generically, generalizing from qq dipoles to 7?.7?.-dipoles (where TZ refers to an arbitrary 

~ ab TZ nP-n J?P TZ t 

representation TZ) one finds their evolution equation to contain the operator [U z \ tr(t U x t U y ) 
(see [51]). The evolution equations of this operator are mapped onto those of 7?.7?.-dipoles or gg- 
dipoles in the two coincidence limits according to 



lim [U z ] ab tr(t U x t U y ) = Cu^tr [U z Ulj , (2.27a) 

TZ TZ a TZ k'kI tz TZ Tit 



lim [U z } ah tr(t Uj U ) = C n tv{U x U y ) . (2.27b) 

2— >y or i y 

The more [/-fields involved, the more constraints are imposed by coincidence limits and one may 
construct whole towers of operators linked downwards by coincidence limits. All these structures 
and relationships are automatically preserved and maintained in full JIMWLK evolution which, at 
least at one loop accuracy, can be simulated numerically. 4 

The theoretical picture can be simplified and the numerical effort required to solve the evolution 
equation can be reduced significantly by truncating the hierarchies. Any such truncation comes at 



4 Limitations arc imposed only by available computational resources with limits on evolution ranges and initial 
conditions that can be accommodated without losing numerical accuracy. 
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the price of introducing an additional approximation. The most widely used truncation of JIMWLK 
evolution is known as the BK approximation. It assumes the factorization 



(S xz S zy )(Y) -+ (S XZ )(Y) (S zy )(Y) , (2.28) 

which turns Eq. (2.25) into a closed equation in terms of (S xy )(Y) only and thus decouples the 
rest of the Balitsky hierarchy. The BK truncation is valid and is parametrically justified in the 
large- N c limit for scattering on a large dilute nuclear target. Using (2.28) in (2.25) we obtain the 
BK evolution equation 

d :(S xy )(Y) = ^ [d 2 z K xzy \{S XZ ){Y) (S zy )(Y) (S xy )(Y)} . (2.29) 



dY" " y/v ' 2ir 2 

An alternative truncation has been discussed in [27,51,60] and dubbed the Gaussian truncation 
(GT) in [ ]. In spirit, it approximates the JIMWLK average with Glauber-iterated two gluon 
t-channel exchange with the target. 

This leads to explicit expressions for multi-L^-correlators, in terms of a two point function Qy. X y, 
for example 

(tr(U x U y ))(Y) = d n e - c - e -- , (2.30a) 

{[U z ] tr(i U x t U y )){Y) =C n d n e~~ (ev,-.+ev,.„-ev,. B )-c w ev,-„ (2.30b) 

and a consistent description of evolution for dipoles in all those arbitrary representations in terms 
of a single equation for Q (see [51]) 

d -,GY, xy = - 9 UzK xzy (x-e-^- )) , (2.31) 



dY* Y ' XV ~ 7T 2 J XZV y 

irrespective of the representation 1Z. One may think of GT as a truncation that, compared to 
BK, includes the minimal subset of 1/N C suppressed contributions needed to restore the group 
theoretical coincidence limits (2.27) which are automatically satisfied by full JIMWLK evolution. 
It can be shown [51] that the dynamical content of GT and BK are in fact the same in the sense that 
replacing qq- and gg-dipoles appearing in Eqs. (2.30) (directly or in certain limits) by their large- N c 
counterparts (proportional to e^~^^ Y ' a>v and e~ JVc ^ y> » respectively) maps the BK equation (2.29) 
onto (2.31) and vice versa. GT improves over BK not in terms of dynamical content, but in the 
way this content is mapped onto different correlators of the theory. 

It was shown in [51] at one loop accuracy, that the modification of the truncation encoded in 
Eqs. (2.30) and (2.31) under the name Gaussian truncation leads to slightly better agreement with 
full JIMWLK evolution than the BK truncation. We will use both BK and GT at NLO accuracy 
to compare to data below and will see that GT results in a slight improvement of the fit in keeping 
with the slightly better match of GT with JIMWLK evolution. The main advantage of the Gaussian 
truncation is that it allows to consistently describe general n-point functions such as that on the 
left of (2.2) before the local limit is taken. Unlike the BK approximation it is versatile enough to 
allow us to test the reliability of the phase space approximations that are built into (2.9) without 
being hampered by 0(1/N%) corrections. 
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2.3.2 Evolution at NLO 



At NLO the picture gets even more complicated. At this accuracy not only single gluon Wilson 
lines are added to the original dipoles at the level of Eq. (2.25), also insertions of nonlocal operators 
such as ti (t a U Z2 t b Ul 2 ) (in the case of quark contributions) appear on the right hand side. It is 
clear that a full JIMWLK treatment becomes more and more costly and suitable truncations more 
and more of a necessity. In the light of the increasingly complicated insertions one would expect 
a strict leading l/N c BK approximation to become rather crude. The Gaussian truncation on the 
other hand should remain a viable candidate for a useful truncation. Here we choose an NLO 
treatment that, for the total cross section (2.9) and the diffractive cross section (2.18) allows us 
to use both BK and GT with comparable accuracy, and only keep in mind that once accuracy is 
high enough to consider testing the validity of (2.2) for the total cross section or more differential 
observables as discussed in [50], the Gaussian truncation becomes the tool to choose. Our main 
reason not to use the full dipole evolution as presented in [57] is numerical efficiency. While it is easy 
to include running coupling corrections according to [54-56], an implementation of the conformal 
corrections is impractical. Instead we adopt to substitute the conformal corrections with an energy 
(or rather longitudinal momentum) conservation correction as suggested by Gotsman, Levin, Maor 
and Naftali [61]. This is an attempt to resum DGLAP type corrections that enter small- x evolution 
at NLO that resum collinear contributions to all orders but should not lead to any double counting- 
conflicts with the resummation of running coupling corrections. 

In the BK-truncation the NLO equation to solve takes the form 

^jy^ Y '- x y = 2^2 J d 2 Z M xzy {l — -^^{SY-xzSy-zy ~ Sy-xy) (2.32) 

while its GT counterpart reads 

jLe-CfSy-.w _ fd 2 zM (l- — ^1 - e -^(5v,»«+5r,^-5v,«v) V-^5y„ 

dY - 2 ^J azAnxzy V dYJ\ J 6 

(2.33) 

The kernel function M xzy = K-xzyR^zy is a product of the leading order BFKL/BK kernel lC xzy 
of (2.26) and what one may call the effective running strong coupling R c xzy . The energy conservation 
corrections are represented by the derivative term on the right hand sides. Without it, both 
equations can be solved by a single step of numerical integrations based on the input of Sy or Qy 
alone. To access the derivative one needs to know these functions at two Y values, Y and Y + AY, 
which forces us to use a much more costly iterative procedure described in App.C. 

To arrive at a precise form for R c xzy one should be aware that there exists no canonical way to 
separate running coupling corrections from the conformal contributions at NLO. Only their sum 
is unambiguously defined, to split them apart one is forced to introduce a separation scheme as 
discussed in [54-56,62]. 

For data comparison we will adopt the separation scheme that subsumes most of the known NLO 
corrections into the running coupling contribution as suggested by Balitsky [56] instead of the 
scheme originally suggested in [54,55]. We introduce the shorthand notations 

r=\r\ = \x-y\, r x = \r x \ = \x- z\ , r 2 = \r 2 \ = \y - z\ , (2.34) 
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to refer to coordinate differences 5 and 

M 2 (r) = ^; C 2 = 4e- 5 / 3 - 2 ^ , (2.35) 

to refer to the scales in the coupling constants 6 . Now we may write the running coupling to be used 
in Eq. (2.32) as 

Rf yz =R* s (r,r 1 ,r2) = a s {n{r)) 

which is a combination of three individual running couplings given by the standard perturbative 
one loop running 

4?r 1 

Qs(M)= A^ ln(^/A2) ; A) = (HJVc - 22V»/3 ; N f , N c = 3 => (3 = 9 . (2.37) 

The choice to subsume as large a part of the NLO corrections into the running coupling has a 
large effect if all other NLO corrections are ignored. Comparing evolution the separation schemes 
of [54, 55] and [56] under these conditions has a large effect on the evolution speed as shown in 
Fig. 3, left panel. For completeness this figure also includes the effect of parent dipole running used 
widely in the literature. This procedure postulates an effective coupling R c s := a s (fi(r)) with the 
scale factor C 2 = 4. This choice leads to almost the same evolution speed than Eq. (2.36) as can 
be seen from the leftmost of Figs. 3. 

The idea behind the energy conservation correction is to include an all orders resummation of 
collinear corrections that start to contribute at NLO - they should make up the bulk of contri- 
butions not yet included. Unfortunately, we lack a derivation that would allow us to interlink its 
treatment with the separation schemes used to define the running coupling corrections: double 
counting of contributions can not be excluded without this information. Correspondingly, the en- 
ergy conservation corrections do not alleviate the difference between evolution speeds induced by 
the different running coupling schemes once included in the calculation. This is shown in Fig. 3, 
middle panel and should be considered a major uncertainty in the present approach. However, this 
is of little consequence for data fits. Once we allow A to become a fit parameter, the difference 
of these treatments can be reabsorbed by a rescaling of that factor: after rescaling, the shapes of 
dipole amplitudes differ very little and allow fits of equal quality. We will find, however, that the 
energy conservation correction is a prerequisite to obtain a good fit in the pseudo-scaling region. A 
treatment that explicitly removes double counting from the outset would be much preferable. 

One unavoidable complication remains in the region where R S (Y)A is near one, the region where 
the IR safety arguments that render J1MWLK evolution a self consistent procedure do not apply. 
JIMWLK evolution is justified where R S (Y)A -C 1, since at distances larger than the correlation 
length R S (Y), the correlator part of the evolution equations (2.32) or (2.33) very quickly approaches 
zero, so that for small enough R s any sensitivity on how one regulates the Landau pole in the running 
coupling modified kernel M. xzy disappears - IR uncertainties are effectively eliminated. 

5 r is the size of the parent dipole (qq dipole) and n,ra refer to daughter dipoles (qg and qg dipoles). All 
expressions encountered are fully symmetric in interchange r\ <-> ri. 

6 7b = 0.5772... is the Euler-Mascheroni constant. The scale factor C 2 = 4e -27E-5 / 3 is specific to the MS 
scheme, it would be replaced by 4e -27E in the V-scheme. 



r{ ( a s (fi(r 2 )) 



r 2 \a s (jjt(ri)) 



ri ( a s (/z(n)) 



r 2 la s (/i(r 2 )) 



(2.36) 
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If, however, R S (Y)A is near one, regulator effects become visible. They affect mostly the large 
r part of the dipolc amplitude and integrated quantities such as the evolution speed, as defined 
in (2.39). 

Since, as discussed below, the data lie in a range where the scale separation between A and Q s is 
not safely established, we have to face uncertainties induced by our choice of regulator, our model 
for the IR behavior of the QCD coupling. The uncertainties encountered are shown in the panel 
on the right of Fig. 3 and involve an APT regulator [63, 64] and the extreme case of a practically 
unregulated coupling (corresponding to a cutoff treatment where a s is the running coupling frozen 
only after it reaches 30.). This serves to illustrate the potential size and range of influence of the 
uncertainty induced by the presence of the Landau pole. 




R A R A r = I x-y I | GeV ' ] 



Fig. 3: Left: The evolution speed A as a function of dimensionless units R S A for two different approaches 
mentioned in text. Middle: Evolution speed after adding in the energy conservation correction on top of running 
coupling Right: The running couplings (regulated vs. unregulated) used in different BK schemes plotted against 
the parent dipole size r. The vertical lines bracket approximately the region of saturation corresponding to HERA, 

0.8 < Rs < 3.7 for 5 x 1(T 7 < x bj < 0.02. 



2.3.3 Scaling behavior and self-consistency at NLO 

JIMWLK evolution as well as both its BK and GT truncations are IR safe in the sense that in a 
transversely infinite medium and an initial condition with short enough initial correlation length 
R S (Y ) all further evolution is governed by contributions on perturbative scales with most contri- 
butions arising near R S (Y) < R s (Yq). Under these conditions the integrand vanishes exponentially 
near the Landau pole and one may argue on physics grounds that contributions to evolution at the 
Landau pole may be neglected. In practice this is achieved by introducing some form of regulator. 
While at asymptotically high energies where Q s ^> 1 GeV the precise choice of regulator cannot 
affect evolution, real data are far from that region and the choice of regulator will impact any 
comparison with experiment. For definiteness we choose what is known as an APT regulator, see 
for example [63, 64] 7 , and simply subtract the Landau pole from the expression for the running 
coupling. The advantage of using an APT regulator over the other possibilities is that it produces 

7 APT stands for Analytic Perturbation Theory. We use only its treatment of the Landau pole a practical way to 
regulate the running coupling. 
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a very smooth behavior of the effective running coupling i? ff even if it depends on all three scales 
jtt r , /x ri and fj, r2 . If one were to simply freeze the individual couplings below some fixed scale in 
i? e ff of Eq. (2.36) the result would be discontinuous. This is partially compensated by a small JC 
in the evolution equation as long as energy conservation corrections are ignored. The stability of 
the iteration procedure used to implement the energy conservation correction described in App C, 
however, is strongly reduced by such discontinuities. The restriction to transversely infinite media 
is conceptually more serious: The JIMWLK framework does not correctly describe the transverse 
growth of a finite target. Near the edge, gluon densities become small and the evolution equations 
match up with the BFKL equation, unphysical Coulomb tails are no longer shielded by a density 
induced correlation length and nonperturbative contribution start to dominate. For this reason one 
usually uses JIMWLK equations to describe the in medium behavior and models any size effects. 
The simplest possible treatment for dipole cross section would be to postulate a separation of r- 
and 6-dependence in simple correlators like the dipole amplitude of Eq. (2.5) and introduce a profile 
function T(b) to model the dipole cross section as 



a qq (Y, (x - yf) = „ Nf. xy with Nf. xy = N% m _ y and a := 2 / d 2 b T(b) . (2.38) 



With this assumption, all the nonperturbative information is encoded in a single constant ao, which 
is then matched to data. While such a treatment is adequate for a LO treatment of the total cross 
section, we will see below that such an approach induces uncertainties already at leading order once 
more sophisticated observables are considered, rapidity gap events and the cross sections of [50] are 
among those. See Sec. 4.2 and the discussion leading up to it. 

The numerical impact of NLO corrections is best illustrated by studying its impact on evolution 
speed A, defined as 



This definition is equivalent to the "naive" definition of A as the rate of change of the saturation 
scale, X(Y) = \n(Ql(Y)), wherever strict scaling holds [65], but provides a useful generalization 
wherever strict scaling does not hold: This includes the pseudo-scaling case encountered as one 
steps beyond LO [58]. Where it becomes necessary to plot energy dependent quantities we often 
choose to plot against R S (Y) as the intrinsic ^-dependent scale of our simulations. For simplicity 
we take R S {Y) from 



Fig. 4 shows the behavior of the evolution speed as one incorporates NLO corrections: At leading 
order, i.e. with fixed coupling and without the energy conservation correction (Fig. 4, left) different 
initial conditions in the course of evolution towards higher Y (moving along the curves from right to 
left in the figures, from large R e (Y) to small R S (Y)) all merge up with the asymptotic scaling regime 
in which A becomes a constant. As running coupling is turned on (Fig. 4, middle), true scaling turns 
into pseudo-scaling: even the asymptotic speed remains V-dependent. For R S (Y)A < 1, evolution is 
slowed down drastically. Adding our last NLO ingredient (the energy conservation correction) leads 
to further slowdown as shown in Fig. 4 on the right. For comparison, all the panels in Fig. 4 show 
the A range that leads to successful fits in our own fits and the scaling models of GB-W and IIM. 




(2.39) 



S q y q xy (\r\ = Rs(Y)) 



1 



(2.40) 



2 



17 




10"' ltf 2 10"' 10° 10 1 10" 3 10" 2 10 1 10 D 10 1 10 3 10" 2 10"' 10° 10 1 

RA RA RA 



Fig. 4: Left: LO i.e. fixed coupling BK evolution: asymptotically the evolution speed A — > const x a 3 for 
any relevant initial state. Middle: NLO BK evolution: the running coupling slows down the evolution. After 
the initial state effects are erased, the evolution speed settles on the same asymptotic line. Right: DGLAP type 
NLO corrections slow down the evolution further. 



One concludes that NLO corrections to evolution (running coupling contributions in particular) are 
an essential ingredient to a successful fit to HERA data. 

In preparation to the discussion of the fit procedure below let us note that we have intentionally 
plotted evolution speeds against R S (Y)A, where, at this point, A is the QCD scale in the coupling. 
In the fits below, however, we will allow A to vary. The reason for this is twofold: we have observed 
already that such a change of scale can absorb most of the differences between different coupling 
separation schemes, beyond that one may also think of this a resummation of nonperturbative 
effects in the initial condition that affects both the values of R S (Y) at all Y and the associated 
evolution speeds X(Y) (through the size of the coupling). If the energy conservation correction is 
omitted as in [45], the latter is main justification to treat the scale factor as a fit parameter. In any 
case, it is through this fit procedure that we set the overall scale of R S (Y) in physical units. 8 

Evolution with NLO effects included allows us an additional cross check on the self consistency of 
our tools: For the calculation to be self consistent, there should be a clear hierarchy of size between 
leading and subleading contributions, unless the subleading contributions introduce a qualitatively 
new feature. This does apply to our discussion of evolution speed near the scaling regime: a 
new qualitative feature (scale breaking and the appearance of running coupling effects) change 
evolution speeds dramatically, while the effect of energy conservation corrections induce only minor 
corrections. This, however, is not universal: the relative size of the energy conservation corrections 
compared to evolution without it strongly depends on the shape of the solutions. Near or in the 
pseudo-scaling regime the energy conservation corrections have a relatively small effect. Far from 
the pseudo scaling regime, the energy conservation correction dominates the r.h.s of the evolution 
equation over a large range of R S (Y) scales. This difference in behavior is shown in Fig. 5. It is 
accompanied by severe numerical stability problems away from the pseudo-scaling region: We have 
employed iterative methods to obtain the derivative term as well as backwards difference methods 
(once two or more adjacent time steps are known with sufficient accuracy) and both require step 
sizes Ay beyond anything even remotely practical to stabilize the numerical results 9 . 

8 [45] fix the value of A and introduce a scale factor c to achieve the same goal. 

9 This is also the reason why we have chosen an approximate iteration procedure outlined in App. C to construct 
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Fig. 5: Relative size of the energy conservation correction: the contributions to the r.h.s. ofEq. (2.33) at a fixed 
R S (Y) (left: l/i? s <C A, right: i? s near A) are split up into the contribution without the energy conservation 
correction, labeled [R.H.S.] r c (it contains only the running coupling corrections) and the energy conservation 
correction -^[R.H.S.Jrc. In the scaling regime, the the energy conservation correction is subleading. Away 
from the scaling regime (exemplified by a Gaussian correlator shape at the same correlation length) the energy 
conservation correction dominates. 



Both these observations, the dominance of the energy conservation correction and the numerical 
stability, lead us to believe that the evolution equation in its present form is less reliable far away 
from the scaling regime. It would appear that additional resummations are necessary to reliably 
address the region far from scaling. The nature of such corrections is not known at present. We 
therefore advocate the use of the pseudo-scaling regime in data comparisons as long as fits in this 
region are at all possible. 



3 Lessons from the total cross section 
3.1 General features 

In the following we will first confront the data with evolution assuming three light quarks with 
current quark masses < 5 MeV. This is sufficient to discuss the main features of the fit and any 
fit tensions. Fits with only a single mass are much faster to do and hence a more efficient tool to 
clarify such systematic questions. 

After scouting the terrain in this manner we verify that the inclusion of quark masses will not 
change our conclusions and explore which phase space ranges are most affected by the inclusion of 
mass effects. 

One generic feature is shared by all fits we have undertaken: the values of R s appear to be well 
constrained by data, an impression of this is shown in Fig. 6 which shows alternatively Q S (Y) = 

the solution in the pseudo-scaling region. 
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1/R S (Y) and Q S (Y) — 4/R s (Y) as they emerge from the fits overlaid on the data used. Sloped 
dotted lines indicate where the qq-dipole amplitude crosses 0.1 and 0.9 respectively. Details on these 
fits will be given below. With 1/R S (Y) of the order of 1 GeV one would expect that nonperturbative 
contributions to evolution are inevitable, but to judge nonperturbative influence in a meaningful 
way, we need a more detailed analysis. Below we will comment on two aspects: IR effects in the 
evolution and phase space features in the case of mass effects. 




10" 2 10"' 10° io' 10 2 10 2 10 4 10" 2 10"' 10° io' io 2 io 2 io 4 
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Fig. 6: Left: The phase phase of all (inclusive and diffractive) HERA data included in fits together with 
the saturation momenta Q 2 (Y) = l/Rg(Y) obtained from different fit approaches. The horizontal lines at 
x — 5 x 10~ r and x = 0.02 indicate the small x range used in the fits. In diffractive fits no restrictions were 
used. Right: To illustrate that slight changes of scale definitions strongly affect the appearance of this plot we 
replace Q 2 S {Y) = 1/R 2 S (Y) by Q 2 S (Y) = 4/R 2 (Y). 

The second feature that is quite well constrained by data, irrespective of details of the theoretical 
input, is the slope of R S (Y) in Fig. 6, or alternatively the evolution speed A(Y). Experience from 
GB-W and IIM model fits with A = .31 and A = .29 respectively as well as our fit experience 
with solutions from the evolution equations all establish a viable range for X(Y) that falls into 
the successful fit range indicated in Fig. 4. If this is generic, fits without energy conservation 
corrections will necessarily have smaller R S A values than fits that include such a contribution to 
evolution (Fig. 4, right). Since the physical values of R S (Y) are very well constrained this implies 
that fits without the energy conservation correction require a noticeably smaller value for A that 
fits that include energy conservation corrections. The former get the contributions from a region 
with smaller coupling and are less sensitive on the choice of IR regulator for the Landau pole. 
The considerably smaller A values necessary without the energy conservation corrections, however, 
would indicate a larger nonperturbative resummation entering the evolution equation: none of the 
alternatives is free of nonperturbative effects. 

In an actual fit with solutions of an evolution equation, one needs to choose an initial condition, 
in our case always of the form of the GB-W model, specify the initial correlation length in units of 
A and evolve this initial condition in Y. As we proceed solving the equation the correlator shapes 
change away from the GB-W form and approach the asymptotic or pseudo-scaling regime. R S (Y) 
shrinks and A traces a trajectory as indicated in Fig. 4. For any given trajectory one then needs to 
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relate the rapidity variable of the simulation to the physical rapidity by determining what we call 
Offset or Y Q g to pinpoint Y — 0. 10 On a given trajectory, Y a g is most closely linked to evolution 
speed. The fit of Y Q g is done simultaneously with a fit of the physical units on R s (Yq = — ln(2-10 -2 )) 
by varying A, and a fit the overall normalization er from Eq. (2.38). 

A fit in the pseudo-scaling region, where the shape of the dipole amplitude is fully determined by 
the nonlinearity with details of the initial condition erased everywhere but in the extreme UV, the 
fit is therefore a three parameter fit in terms of Y Q g, A, and <7n. 

Away from the asymptotic pseudo-scaling region many additional features of the initial condition 
survive and may affect the quality of the fit. At present we have no systematic tools to scan the 
space of initial conditions and very little constraints from theory. Our efforts below are meant to 
shed some first light on the issue using a very hands on attitude that is limited in scope mostly by 
the cost of creating a single trajectory on the one hand and the search of fit parameters for a given 
trajectory on the other. 

3.2 Systematics from fits with light quarks, energy conservation included 

Before we turn to fit quality, let us first collect the main features the data fits impose on correlator 
shapes and evolution speed in Fig. 7. 




Fig. 7: Left: The evolution of Ny q (r) in different approximations and models. The data align JIMWLK 
based descriptions and the MM model in the region N^?{r) < 0.4 - the GB-W model deviates significantly. 
Middle: The asymptotic (perturbatively imposed) solutions of the evolution equations are steeper after the 
energy conservation corrections are included. This benefits the fit quality overall and at large Q 2 in particular. 
Right: The evolution speed A for all JIMWLK approximations considered as a function of correlation length R s . 
Error bands are extracted from the condition x 2 /dof = 1, see Appendix B. 



We have already stressed at the outset of Sec. 3, the most tightly constrained feature in all the fits 
is R S {Y) in physical units. This also manifests itself in the fact that after the fit is performed, the 
shapes of dipole amplitudes agree very well within an order of magnitude below R S (Y) no matter 

10 Since we only start compare to data for x < 2 ■ 10~ 2 we may well obtain negative values of Y ff and still have 
dipole cross sections to cover all the data range considered. We have no ambitions to extend or parametrizations to 
larger x. 
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if the calculation is based on an evolution equation or a model. This is shown in Fig. 7, on the 
left panel for BK, GT, BK with parent dipole running (labeled BK p( j) and the two models GB-W 
and IIM. The models show noticeable deviations from the fits based on the pseudo-scaling solutions 
from evolution equations at |r| > R S (Y). 

We had already seen in Fig. 4 that the main effect of the energy conservation correction is a further 
slowdown of evolution above that induced by the running of the coupling. The middle panel in 
Fig. 7 illustrates the second major impact by comparing BK and GT fits (which include the energy 
conservation correction) and an asymptotic BK fit without this correction labeled [BK] rc : correlator 
shapes in the asymptotic regime stay steeper than in the corresponding evolution with the energy 
conservation correction omitted. 

Evolution speeds corresponding to the fits are shown in the rightmost panel in Fig. 7. We indicate 
an associated error band only for the GT fit. It encloses all fits that include the energy conservation 
correction, but is clearly separated from the fit without the energy conservation corrections [BK] rc , 
complementing the shape deviation already observed. The two models, by construction, have 
constant evolution speeds that fall near the average of evolution speeds obtained from the BK and 
GT fits, confirming our expectations. 

Table 1 provides an assessment of the quality of the fits illustrated previously in terms of dipole 
amplitudes and evolution speeds in Fig. 7. The table reflects our fit strategy of first obtaining a 
fit at low Q 2 < 45 GeV 2 which we then attempt to extend to a larger Q 2 range up to 1200 GeV , 
monitoring any change in fit parameters required in the process. In the larger Q 2 range the models 
(GB-W and IIM) are known to fail, they were designed only with small Q 2 values in mind, but also 
the small x approximation employed in the evolution equations will have to break down eventually 
- even with partial resummations like the energy conservation correction built in. The success of 
such an extension is a measure as to how efficiently the resummations recapture large Q 2 effects. 

x bi < 0.02 BK pd BK GT GB-W IIM [BK] rc 

A(a; bj ) 0.27-0.34 0.26-0.35 0.26-0.35 0.31 0.29 0.31-0.44 

*bj 6 [jjjr-,0.02]) 

A [MeV] 

Q 2 < X 2 / 224 
45 GeV 2 CT [GeV 



-2i 





98.7tg;I 


104.7±f?;| 




2-0 — iq4 


50.4 


0.818 


0.811 


0.810 


1.401 


0.828 


1.760 


54.01 


55.05 


55.33 


44.59 


51.47 


56.84 


0.979 


0.812 


0.805 


1.978 


1.037 


4.783 


53.69 


55.01 


55.35 


44.48 


51.08 


54.79 



Q 2 < X 2 / 2 95 
1200 GeV 2 <T [GeV -2 ] 

Table 1: Fit results to inclusive data from [34-36,41]. x 2 /dof is below one for a wide range of A values 
indicated by the errors listed, see Fig. 19. 



The best fits are obtained when we use the NLO evolution equations once energy conservation is 
included (columns labeled BK and GT as well as BKpd): 11 The fit quality is excellent over both 
the small and the large Q 2 range. Already the low Q 2 range (Q 2 < 45 GeV 2 ) determines both the 
physical units of i? s (V) and the ideal evolution speed A(Y) (via Y g). To extend the fit to the full 
Q 2 range covered by the data below x = 2 • 10~ 2 only cr needs to be readjusted. Note the excellent 

J1 A11 these fits are performed with light quarks only, the role of physical quark masses is discussed in Sec. 3.3. 
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X 2 /dof values over the whole Q 2 -range for both BK and GT truncations in the pseudo-scaling 
region. 

An asymptotic fit with the energy conservation omitted is clearly unworkable: its x 2 -value is barely 
acceptable already for Q 2 < 45 GeV 2 and indicates an outright failure in the broader Q 2 region 
up to 1200 GeV 2 - it fares worse than the models (if extrapolated into this region) by far. The 
origin for this is the shallower shape of the dipole amplitudes observed in Fig. 7 (middle panel). 
The data require steeper correlators and lower evolution speeds. Evolution speeds do slow down 
as evolution proceeds to smaller R S (Y)A, however, at the same time correlators flatten out further 
to asymptotically approach the shapes of the fixed coupling case as the running of the coupling 
slows down with shrinking R S (Y)A values. This leaves a very small window for A (viewed as a fit 
parameter) in which the correlators are still tolerably steep, but evolution speed is already small 
enough. As a result all features of this fit deviate from those shown in the left panel of Fig. 7: 
evolution speed is still quite large, and the match of the dipole correlators below R s observed in 
the left panel of Fig. 7 is lost as well. 

The success of the asymptotic fit with energy conservation included could be interpreted as an 
indication that the inclusion of the energy conservation correction gives a better match to the 
perturbative anomalous dimensions that govern the large <3 2 -behavior. This is in contrast to the 
preasymptotic fit of [45, 4G] (which omits energy conservation corrections) where relics of the steep- 
ness of the initial condition -non-perturbative in nature- allow for a good fit quality at all Q 2 . We 
will provide more details on this comparison in Sec. 5. 

The models, GB-W and IIM, are clearly limited to a smaller Q 2 region. This holds even for the IIM 
model, which does incorporate additional perturbative information in the form of BFKL anomalous 
dimensions. 

Fig. 8 shows asymptotic solutions and models against a subset of data to illustrate fit quality 
in different Q 2 ranges. The asymptotic fits with energy conservation included remain valid to 
astonishingly large Q 2 values, exceeding the 1200 GeV 2 range over which we have kept track of 
X 2 -values above. This generic picture is reinforced if we contrast the scaling behavior of theory 




Fig. 8: F2(x,Q 2 ) as a function of x from a fit with light quarks compared to a subset of data . Different 
descriptions are qualitatively the same at low and moderate values of Q 2 (on the left and middle respectively). 
At high Q 2 (on the right) only the BK and GT fits compare well to the data. 

and data as done in Fig. 9 for GT (left), BK (middle) GB-W fits (right). Shown are j*p cross 
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sections (see Eq. (2.9) for the connection a 1 * p Fa) as a function of scaling variable Q 2 /Q 2 {x^) 
for the three cases after the fit is performed. The borders between moderate and high virtualities 
(above Q 2 = 120 GeV 2 ) are roughly indicated by the arrows overlaid on the plots. Data sets from 
different running periods are plotted by using different colors and symbols 12 . To avoid overlapping, 
the fit results are separated from the experimental data by dividing the normalization factors 00 
out. Note that the saturation scales are different functions of Xbj in each of the panels, although the 
differences between the two truncations (GT, left and BK, middle) are so small that no deviation 
can be discerned visually: both slopes coincide excellently with the experimental data up to the 
largest Q 2 /Q 2 (xbj)- Contrary to that, the GB-W model shown on the right can not resolve the 
high Q 2 data: the slopes of data and theoretical predictions start to deviate at large Q 2 /Q 2 (xbj). 
A more detailed picture of the fit quality is given in Fig. 10. 




,, |( l l ., sfi lU,, |(1 .., sfi lU, Sfi 

10" 10" 10 10 10 10" 10" 10 10" 10 10" 10" 10 10 10 10" 10" 10 10" 10 10" 10" 10 10 10 10" 10 10 10" 10 



<w,<y « !/< 3;<v q 2 *#v 

Fig. 9: The 7*p cross section as a function of scaling variable Q 2 /Q 2 S . The plots show both data and theoretical 
predictions (shifted downwards by 00 to avoid clutter). Left: GT. Middle: BK. Right: The GB-W-model. Data 
from [34-36]. 



3.3 Quark masses 

As already mentioned, quark masses are formally a subleading effect from the perspective of our 
small x resummations, but they do impact final state phase space and the width of wave functions 
in impact factors quite severely. The use of constituent quark masses of 140 MeV for light quarks 
and 1.4 GeV for charm is quite widespread in the context of the GB-W model and typically used 
in a restricted Q 2 -range (below 45 GeV 2 ). Since quark masses by nature are a nonperturbative 
feature that should have its main effect at small Q 2 one should expect that the inclusion of quark 
masses will not spoil the excellent fit quality that was obtained in our light quark fits for Q 2 above 
the heaviest quark included. From diffractive measurements we know that charm quarks should 
contribute significantly to the HERA cross sections while bottom quark contributions are negligible. 
One should therefore complement the three light quarks used above with a charm quark. This brings 
in a quark mass of 1.2 — 1.4 GeV depending on whether one considers current or constituent quarks, 
both of which are of the same order as 1/R S (Y) in physical units and thus one would expect at least 
some complications in the nonperturbative sector. A straightforward fit with three light and one 

12 The same notation for the inclusive data is applied throughout. 
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heavy quark shows that the inclusion of the charm quark reduces fit quality mainly in the Q 2 range 
below 1 GeV 2 (First two columns in Table 2). This is about as far as we can go in our analysis 
without introducing any model elements that modify the low Q 2 behavior in some ad hoc manner. 



Shaded band, Y: 3.91-8.52 




Fig. 11: Cross sections including charm quarks with m u ^,a,c = {3, 5, 105, 1270} MeV for Y = {5, 10, 20, 40}. 
The curves are calculated with the parameters corresponding to Xeft = x(l + Ql/Q 2 ) (or equivalently Y e ff = 
Y — In (l + Ql/Q 2 )) in Table 2. Left: o%S p and rf* p (blue dashed) plotted separately. Below the largest quark 
mass used, a starts to drop. This is incompatible with current conservation in the limit Q 2 — > and reduces 
the quality of the fit. Right: The charm fraction (j^* p /(T* * p as a function of Q 2 . The horizontal line indicates 
the large Q 2 limit e 2 /^2j.e 2 — 2/5. The data are from [66-68]. Agreement is clearly qualitative at best. 



We should, however, at least qualitatively discuss a nonperturbative modification introduced by 
Golec-Biernat and Wiisthoff [30] , a modeling device to accommodate nonperturbative contributions 
at small Q 2 . They have suggested to evaluate the dipole cross section at 

XcS ~ X Q2 - Q2 + W 2 V- 1 ) 

instead of x in order to guarantee VF-independent cross sections at small Q 2 as required by current 
conservation in the photo-production limit at Q 2 = 0. We have found that with the data set used 
here, (3.1) in fact improves the fit for Q 2 < 45 GeV 2 in the GB-W model. When using asymptotic 
solutions to the evolution equation, however, fit quality goes down (taken over the full Q 2 range) 

o 2 

as shown in the third column of Table 2. In fact even the "inverse" modification x c fr = x — r 

works better from a % 2 perspective (fourth column), despite a clear lack of supportive arguments. 
One possible reason for the failure of (3.1) is conceptual: it is not compatible with the factorization 
into impact factors and Wilson line correlators that is at the core of our renormalization group 
picture. Quark masses are properties of projectile constituents and have every reason to show up 
in impact factors. On the other hand, it is hard to imagine that they should be resummed into a 
feature of the energy dependence of the Wilson line correlators, which are purely determined by 
target properties. But this is exactly what is done when using (3.1). From this perspective the 
only scale available for use in a modification at small Q 2 is in fact Q 2 {x). We tentatively suggest 
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to replace (3.1) with 



Q 2 + Ql(x) Q 2 + Ql(x) 

X ^ = X —Q, = Q2 + W 2 ■ (3-2) 

While this does not freeze the cross sections at fixed W, it does flatten them noticeably at small Q 2 
and leads to an improvement of the fit quality, both with a schematic mass pattern of m u ,d,s,c — 
{5, 5, 5, 1400} MeV and the current quark mass pattern m Ut d,s,c — {3, 5, 105, 1270} MeV as shown 
in the last two columns of Table 2. We emphasize that (3.2) is only a conjecture to resum part 
of the nonperturbative contributions below Q s , one should not take too much encouragement from 
the improvement of \ 2 alone, see the ad hoc success of the redefinition in column 4. 

Q 2 +4m 2 f 



XcS 


.r 


X 

(Q 2 > 1 GoV 2 ) 


x Q2 


X Q 2 +4m 2 

ad hoc mod. 


x Q2 


x Q2 

(phys. masses) 


V off 


-5.86 


-6.87 


-5.63 


-5.99 


-6.50 


-6.73 


A [MeV] 
a GeV~ 2 


88.1 


63.5 


89.7 


90.5 


70.3 


66.1 


57.48 


73.04 


61.20 


51.69 


68.94 


70.52 


X 2 /dof 


1.41 


1.05 


1.59 


1.23 


1.13 


1.14 



Table 2: Including a charm quark; nonperturbative modifications. For the second column data with Q 2 < 
1 GeV 2 are removed, reducing the d of from 295 to 221. Quark masses are m Uy d,s,c = {5,5,5, 1400} MeV in all 
but the rightmost column where m Ut d,s,c ~ {3,5, 105, 1270} MeV. 



4 Lessons from the diffractive data 

4.1 The need for NLO contributions to the impact factors 

Diffractive HERA data extend down to ft ~ .04. Thus overall Y = ln(l/a;) and Y gap = ln(l/xp) 
remain comparable while Vf rag = ln(l//3) remains too small for multiple gluon emission to build up 
within the projectile fragmentation region - contributions from the ggg-component of the impact 
factor (which has its first contribution at NLO), however start to play a role even at such moderate 
ft values as already observed in the pioneering papers of [69,70] and reiterated in [71]. The main 
reason for that is that the qq contributions given in Eq. (2.18) strictly vanish at ft — ¥ 0, even 
with NLO effects to the evolution of the Wilson line correlators taken into account (see Figs. 13). 
Any gluon component in the impact factor, on the other hand, will generate a nonvanishing cross 
section in this region of phase space: the NLO contributions to the impact factors are the leading 
contribution at small ft. 

Unfortunately, no full expression for the c/c/g-component is available, only the large Q 2 and small ft 
limits (without /3-evolution) are known exactly. However, an interpolating form has been suggested 
in [53]. Below we find that nonperturbative contributions related to the target profiles by far 
dominate the uncertainties, and, for simplicity, we content ourselves with the large Q 2 expressions 
of [69, 70] to estimate the contributions. 



27 



The starting point then is 



dp 



7* A— s-Jfp 

dp 



dp 



da 



7* A— s-Xp 

qqg,T 



dd 



(4.1) 



LL(Q 2 ) 



where the corresponding structure functions xpF®^ , xpF®^ , and xrF^J^ [all functions of 

(xp, Q 2 , /?)] are obtained by dividing out a factor 4 ^"^ m . The first two terms in (4.1) are given in 
Eq. (2.18), they contain the LO impact factors with only a qq Fock component in the final state. 
The last term is the large Q 2 part of the contribution of the NLO impact factor: only the transverse 
part contributes, the longitudinal part being of higher twist. At large Q 2 the qqg-Fock state appears 
in a configuration in which the inter-quark-distance is tiny compared to 1/Q S and one may take 
the corresponding coincidence limit. The qq part is then indistinguishable from a gluon and the 
analytic expression may be cast in terms of gluon dipole amplitudes Ny '(r, b). The corresponding 
expression was first given by Golec-Biernat and Wiisthoff [ ]. In it, masses are set to zero, in the 
spirit of a large Q 2 expansion. As with the first two terms of Eq. (4.1), which were already given 
in Eq. (2.18), we present this result in a notation inspired by [53], but with the 6-integral not yet 
performed. This allows us to assess the quantitative impact of our lack of precise knowledge of the 
6-dependence in all three contributions. We retain the assumption that the dipole amplitudes only 
depend on |r| and are independent of the orientation of the dipole. Then 



^T XP 



dp 



LL(Q 2 ) 



a cm a s C f N c f 1 dz 

{l-zf 



tir 2 Q 2 



1 - 



(l-z)Q 2 



dk 2 In 



OO 

x J dr 2 dr' 2 <j) gg (z,\k\,\rl\r'\) J d 2 b (r,b) N°? (r' ,b) 



(4.2a) 



where 



J gg 



:k 4 J 2 (\k\\r\)K 2 



1 - z 



-h 2 



Kn 



1 



ft 2 



J 2 (\k\\r'\). (4.2b) 



Any inclusion of quark masses (or a dependence of the dipole orientation) in such an expression 
would amount to a resummation of subleading effects with little control over their relevance. Sim- 
ilarly, we have no reliable argument to set the scale in the strong coupling a s that appears in the 
prefactor of this expression. Both these assessments are reinforced once one starts analyzing the 
nonperturbative uncertainties in Eq. (4.1), even after it has been updated on the perturbative level 
with full NLO ingredients. To be consistent in our treatment below we will therefore also only 
consider the massless limit for the qq contributions. 

In fact, the largest uncertainty in (4.2) and in the corresponding expressions for the quarks, 
Eq. (2.18), is related to the impact parameter integral which is not under perturbative control. 
Even after using data to set the overall normalization of the total cross section (the main non- 
perturbative parameter entering the LO total cross section), already the leading order diffractive 
contributions (the qq-terms in (4.1)) require additional nonperturbative input. This only gets more 
pronounced at NLO: higher order Fock components in the projectile wave function couple to higher 
n-point functions of Wilson lines, each of which is affected in its own way by non-perturbative 
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effects. This affects the relative normalizations of the terms in Eq. (4.1) as well as the relative 
normalization of total and diffractive cross sections. In practice, this manifests itself in a strong 
model dependence of the normalization of individual cross sections. The main issues here are the 
relative normalization of total and diffractive cross sections on the one hand and the weight of 
individual Fock components such as the qq and qqg contributions in Eq. (4.1) on the other. 



4.2 Nonperturbative aspects of b-dependence and profile functions 

The general attitude for the total cross section - to assume a fixed, x-independent target size, 
which at LO implies taking the qq-dipo\e amplitude as the product of a profile function T(b) and 
r-dependent remainder N^?(r) (normalized to one at r — > oo) leaves us with only a single fit 
parameter, the area resulting from the 6-integration. In the diffractive case, the choice of profile 
strongly affects overall and relative normalization of 

d 2 b Nf(r, b)Nf{r', b) (4.3) 
and 

J d 2 b Np*(r, b)Np>{r\ b) (4.4) 

featuring in the formulae above as well as the more general correlators in their full NLO general- 
izations. 

The simplest treatment would associate a factorized profile with each of the amplitudes, identical 
for both quarks and gluons according to (again with 1Z labeling the representation) 

N$(r, b) ->■ T{b)N?( r ) (4.5) 

and already in this case, the relative weight of diffractive and total cross sections are highly model 
dependent: A box profile Tb ox (^) of height one, normalized in width to produce a factor 

o-o = 2 J d 2 b T(b) (4.6) 

for the total cross section, results in a factor J d 2 bT 2 ox (b) — <tq/2 in the diffractive case. A Gaussian 
profile, which has some phenomenological justification at large |b|, 13 produces an additional factor 
| in the diffractive case compared to the box profile: the area under 2Q auss (6) is half the area 
under TG auss (6). Clearly the relative normalization of the total and the diffractive contributions 
is strongly dependent on the shape of the profile: Arbitrary factors of this sort can already be 
obtained by varying the width and the height of the box profile while keeping (4.6) fixed. While 
one may dismiss box profiles as unphysical, the issue remains: there is by no means a canonically 
prescribed physical profile that would outright eliminate such modeling choices. 

Also a Gaussian profile, justified as it may be in some 6-ranges for quarks, leaves intrinsic uncer- 
tainties: Not only is the overall normalization an issue, but also the relative normalization of parton 

13 This is based on successful parametrizations of meson production data via e~ Bd ^ ^| and constrains the 
shape for |6| > .3fm [72]. 
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species (here quarks vs gluons) is affected. To explore this in more detail, we will invoke Casimir 
scaling, which we expect to hold at least at small 6, as a guiding principle. Assuming the same 
profile factor in front of both quark and gluon amplitudes breaks Casimir scaling for the full b- and 
r-dependent amplitude Sy(r, b) = 1 — Ny(r, b) for effectively all b 2 > and hence, is not compat- 
ible with the notion that the energy dependence is dominated by perturbative gluon emission as 
encoded in the evolution equations and their Gaussian truncation. To guarantee this then requires 
that Casimir scaling is at best weakly broken for central collisions. Models with Casimir scaling 
restored completely are readily constructed: One may, for example, exponentiate the profile in the 
spirit of the IPSat and bCGC models according to 

A^(r,fo)^l-e- c - T ^ e -M . (4.7) 

Alternatively, one might start with a factorized (Gaussian) profile for quarks to define Gy{t, b) for 
general use in the Gaussian truncation via T{b)N Y q (r) =: 1 - e -cvSr (r,b) Tmg aUows tQ calculate 
general n-point functions in the Gaussian truncation with 6-dependence and yields Casimir scaling 
for dipoles according to 

Ny(r, b) ->■ 1 - [1 - T(b)Np(r)] %f . (4.8) 

Taking Casimir scaling as a guiding principle modifies both b and r dependence in a C-r dependent 
manner. In particular it leads to sizable changes in relative normalization of quark and gluon contri- 
butions as compared to the height-one box profile, for which all of these definitions are equivalent. 

Plain exponentiation of the profile as in (4.7) enhances the influence from nonperturbative regions 
of phase space as Y increases: with increasing Y the large |r|-growth of Qy(r) will progressively lift 
up any nonvanishing large \b\ tails of the profile function if they exist at all. This in turn leads to 
y-dependent growth of the overall normalization of the dipole cross section (after the 6-integration 
is done). This results in an inconsistent, unphysical interplay of perturbative gluon emission with 
non-perturbative long range physics. One might attempt to regulate the large r-behavior of Qy{ t ) 
to preclude that, but to do this in a defensible way would clearly require non-perturbative input 
completely outside the scope of JIMWLK evolution. 

The ansatz (4.8) on the other hand does not require any such additional input: large \r\ and 
\b\ behavior decouple. Therefore we use this model below to estimate the impact on relative 
normalizations of cross sections. 

With an eye to the total cross section, we first note that the model leads to representation dependent 
normalization for the 6-integrated dipole cross section. Assuming a Gaussian profile, we find that 
the b- integrated profiles acquire a representation-dependent r-dependence we denote Ny" (r) as well 
as a nontrivial large |r|-normalization which we choose to display explicitly 



/ rf» 6 Ny(r, b) := 2«B d H&)N?(r) ^ **B&&) « 2,B d \\ *» 
J Of Uf 11.6 tor gluons 

(4.9) 

where H(z) :— ipo(l — z) + is the harmonic number of z (expressed via the digamma function 
ipo). Diffractive normalizations differ from the normalizations in the total cross section, and in 
addition, the r- and r'-dependence after 6-integration for a generic representation only factorizes 
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approximately 



with 




(4.10a) 



(4.10b) 



Unsurprisingly, normalizations of quark and gluon contributions differ again: the relative diffractive 
normalizations for quarks and gluons are reliably assessed along the diagonal r = r'. where one 
finds 

/ Sb [N^rM? = 2nB d - H&) ^ 2,B d \\ f , (4.11) 

J Cf Cf II. for gluons 

i.e. the gluon contribution is enhanced by a factor of two compared the quarks. This is not the case 
for the ansatz (4.5) or the use of a box profile of height one. We take this as a practical indication 
that the nonperturbative uncertainties are indeed large and that a precision fit at NLO requires 
refined nonperturbative input or independent phenomenological constraints on the individual nor- 
malization of each Fock component. 

On the practical side one finds that factorization is exact for quarks and a good approximation 
for gluons. For quarks N^{r) = Ny?(r) and (2H(^j) — H{^j-)) = |. For gluons one may, for 

simplicity, even use iVy 9 (r) as inspired by (4.9) to approximate 

J d 2 b N*?{r, b)N^{r\ b) « 2nB d C N 1 * 9 (r)N^ 9 (r') (4.12) 
where C = 1, in line with (4.11). 

These may be used in the momentum space formulae originally derived with height one box profiles 
(which factorize trivially) after properly including the normalization factors - the factorization error 
is much smaller than any nonperturbative uncertainty inherent in the choice of the profile model. 

4.3 Diffractive fits to HERA data 

Given the uncertainties arising both from incomplete NLO impact factors (Sec. 4.1) and nonper- 
turbative aspects of the impact parameter dependence (Sec. 4.2) precision fit to diffractive data are 
out of the question. 

The fits presented below are done with this in mind - mainly to assess if one can get a qualitative 
agreement with data, based on the fit parameters Y D g and A obtained from the total cross section. To 
optimize the diffractive fits, an independent normalization Bd is allowed and its relation to B d ot = 

14 For (4.7) factorization of r-dependence is generically a bad approximation and the result depends strongly on 
the IR regularization at large r. 
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ao/Air is ignored in the fit. How closely the fit results match is treated a consistency check. As 
already indicated, the scale of a s in the qqg component of the diffractive cross section in Eq. (4.2a) 
is undetermined. This factor of a s will therefore be treated as a further fit parameter. While Bd 
and a s must be meaningful in the context of HERA physics, we need a means to accommodate 
the uncertainties exposed above. In this sense, we take the precise values of these parameters as 
scenario dependent. 

In the fits to the total cross section, the differences between the different scenarios using the energy 
conservation correction turned out to be small. For this reason, it is sufficient to consider only GT 
results from our simulations. We retain the GB-W model results for comparison. In what follows, 
the overall normalization of the qq components is the same in all cases since for the quarks the 
profile function factorizes trivially. For the gluon component, the relative normalization is scenario 
dependent and we explicitly consider the following two cases: 

• [fact]: 6-dependence completely factorized: the simplest case given in (4.5). The qq and qqg 
components are equally weighted. 

• [sc]: 6-dependence based on Casimir scaling after using the approximation of Eq. (4.10). The 
qqg component is enhanced by a factor of two relative to the ansatz (4.5), its shape is slightly 
modified. (The qq components remain completely factorized.) 



Ny(r) for the two schemes and the corresponding momentum space amplitudes Ny- 9 {q) (see 
App. D.l for definitions), are illustrated in Fig. 12. The difference between the factorized and 




r [GeV 1 ] q [ GcV J 

Fig. 12: Three gluon dipole inputs: globally factorized N Y 3 (r) Eq. (4.5), Casimir scaling based 6-integrated 
(2H g/4 - H g/2 )2N Y 3 {r) Eq. (4.10) and N^ 3 (r) in coord inate space (left panel) and momentum space (see 
App. D.2 for precise definitions) (right panel). 

Casimir scaling schemes are pronounced: The main effect of [sc] is to increase the effective R s after 
6-integration is done compared to the factorized scheme [fact]. 15 As a result the [sc]-gluon stays 
closer to the quark result than the [fact]-gluon, see Fig. 24 in App. D. 

15 We note in passing that the gluon amplitudes Ny 3 (r) (dashed blue) and N Y 3 (r) (solid black) are almost the 
same: a tiny rescaling in r is required in order to (perfectly) match the curves. 
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Since the literature lists data directly for structure functions XpF 2 



D(3) 



Q 2 /B da~i* A ^ Xv 



, the 



actual task is to compute each component of the cross section Eq. (4.1) with the parameters obtained 
from the total cross section and then to optimize B d and a s in the following expression, 



' qg,T I B d = l + t qq,L I B d = l 



■ a ' x r b W9,T I Bd =l, a .=l 



(4.13) 



The numerical burden of this procedure can be vastly reduced by replacing the coordinate space 
expressions for these contributions shown in the text by their momentum space counterparts, see 
App. D. XfF^ an d its components are visualized in Fig. 13: the overall behavior of x&F% is 
analogous in all cases, including the GB-W model. 





Fig. 13: 

fixed Q 2 
Eq. 4.14. 



Behavior of the contributions given in Eq. (D.2). Left: Change in xpF®^ as ip is decreased with 
Middle: Q 2 is increased with fixed a^p. Right: A comparison of xpa 
See the text for more details. 



D ' 3 ' and xtF^ according to 



Np(r): (4.5) GB-W GT N^{r): (4.10) GT 

B d 2.85 4.70 B d 4.70 

a s 0.54 0.44 a s 0.31 

X 2 /552 1.44 1.44 X 2 /552 1.47 

Table 3: The results of the fits to the diffractive data [37-41]. Left: [fact] Fully factorized 6-dependence of 
Eq. (4.5). Right: [sc] 6-dependence based on Casimir scaling according to Eqs. (4.10), (4.11) and (4.12). 

Fig. 14 presents a subset 16 of x-pF 2 (Q 2 ,f3, xr) data from HI [39] together with the theoretical 
predictions. We show three practically indistinguishable scenarios: 

• [fact]: GT with the factorized ansatz (4.5) (denoted "gt, fact.")] 

• [sc]: GT based on the Casimir scaling ansatz of (4.10) (denoted "gt, sc."). 

• GB-W with the factorized ansatz (4.5) (denoted "gbw, fact.") 



•' Although not shown, the data from ZEUS [37,38,40,41] were also included in the fits. 
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Fig. 14: xpF^ 3 (xpQ 2 , (3) as a function of xp. The data are from HI [39]. 
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GT p < 0.5 p > 0.5 

Nf{r): (4.5) x 2 /264 = 1.58 x 2 /288 = 1.30 
JV^(r): (4.10) x 2 /264 = 1.64 x 2 /288 = 1.31 



Table 4: The diffractive data split into two subsets j3 < 0.5 and /3 > 0.5. The data with /3 > 0.5 fit somewhat 
better (a region where qq components dominate). 

The corresponding fit parameters are presented in Tables 3. 

The parameter a s turns out to be strongly scenario dependent, reflecting the differences in the 
shapes, scales and overall normalizations of the gluon amplitudes shown in Fig. 12. For GT, the 
common normalization Bd as well as the fit quality x 2 /dof are of the same order in all cases. The 
values of Bd are systematically smaller than the experimental value Bd = 7.1 GeV 2 reported 
in [38] but, they are consistent with B d ot = a^jAi: obtained from the total cross section. 17 BK 
evolution (not shown) results in fits of the same systematic behavior and fit quality as the GT fits 
presented here. 

Fit quality varies noticeably across phase space: splitting the data into two subsets P < 0.5 and 
P > 0.5 or equivalently Q 2 < M\ and Q 2 > M|-, 18 we observe that the best match with the data is 
obtained at large P where the qq components dominate, see Table 4. This reinforces our statements 
that a better treatment of the qqg amplitude is required. The presently implemented improvements 
in our treatment of the ^-dependence only lead to a tiny improvements in fit quality: the main 
differences between the gluon amplitudes resulting from the two schemes dubbed [fact] and [sc] are 
effectively absorbed into the normalization of the qqg component via parameter a s . 

As a check of consistency, the fits are also performed by using the reduced cross section (as was 
done in [53]) 



where s ep — 318 2 GeV 2 at HERA. We note that for the present diffractive data y < 0.5 so that ay* 

hardly deviates from over the bulk of the data range 19 . As a consequence, one would at best 

expect a tiny improvement in the fit quality in any approach once the reduced cross section is used. 
We find indeed that % 2 values improve almost imperceptibly to % 2 /552 = 1.39 and x 2 /552 = 1.42, 
respectively. Consistent with this, the parameters remain practically unchanged: Bd = 4.75 GeV~ 2 
and a s = 0.43 for (4.5) and B d = 4.75 GeV -2 and a s = 0.31 for (4.10). 

To pinpoint precisely where the differences arise consider the rightmost panel in Fig. 13: this 
compares oy*^ with for fixed xp and Q 2 as a function of p. It should be noted that if 

17 Note that the use of 7b ox (6) of height one induces an overall factor of two in all components relative to TGauss(fr) 
in the framework of (4.5). Since it can only be absorbed by a redefinition of in the diffractive fits, this would 
lead to major inconsistencies. 

18 As with the data for the total cross section, there is kincmatical correlation in the data range: p is small if 
Q 2 <C M%- in this case, x^j = Xfp is also small and thus s 1 » v Fa Q 2 /xbj is large. Since s 7 » p s ep in the 
experiments is finite, small values of p are more likely paired with low Q 2 . 

19 For 322 points out of 512 y < 1/4 implying that the factor in front of in (4.14) is < 0.04. The data with 

l/4<J/<l/2 are most likely associated with large Q 2 . 




y 




(4.14) 
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the diffractive phase space is presented in this way, one must take into account that the kinematic 
limit y < 1 is violated at small /3. Two conditions y — 1 and y = 0.5 marked by the red circles 
and blue squares overlaid on each curve are indicating the smallest possible f3 for fixed xp — 10 -3 , 
Q 2 = {7.5, f5,30,60} GeV 2 and s ep = 3f8 2 GeV 2 (set by /3 min = Q 2 /(x P s ep y)). The former is the 
absolute kinematic limit of HERA whereas the latter approximately marks the lowest values of (3 
of the existing data. 



4.4 Ratios with the total cross section 

The approximate constant ratio of the diffractive to inclusive cross sections was originally observed 
at HERA by ZEUS collaboration [40]. To compute this observable, we use the relation (see [40,41]) 



1 da^ Xp (Q 2 , M x , x P ) ^ (2tt) 5 



2M X dM x Q 2 {Q 2 + M 



x v F? {s) (Q 2 ,M x ,x r ) (4.15) 

x) 



for the diffractive cross section. The ratio of the diffractive cross section to the total cross section 
can then be computed by 

RdiS( , Im^Mx da^ Xp (Q 2 ,M Xl xp)/dM x 

Rtot ( S 7*p) = ~tat fa. Q2) i Q < S ~t'P < S e P ) ( 4 - 16 ) 

where the integration boundaries, i.e. the bins for the diffractive mass M x , are determined by the 
experimental setup (see Table. 5). 

In the calculation of R^ given above, no new free parameters are introduced. Basically, Rf^f is 
parameter free since the normalizations cancel trivially. However, as mentioned earlier, we ignored 
the value of B d ot = ao/Air from the total cross section and allowed a distinct normalization in the 
diffractive fits. For instance, in the case of GT the ratio of the optimal normalizations is found to be 
Bd/B d ot w 1.07, i.e. fairly close to one. In practice this means that all HERA data can be resolved 
by the same normalization B d ot = uq/At: w 4.40 GeV -2 : a single parameter fit to the diffractive 
cross section then results in a s = 0.34 with a reasonable x 2 /552 = 1.52 for the scenario (4.10). The 
region /3 < 0.5 is responsible of the reduction of the overall % 2 /dof: at j3 > 0.5 the fit quality is 
actually slightly improved. 

In Fig. 15 are shown the latest i?tot data from [41] together with the theoretical predictions based 
on the fits to the inclusive and diffractive cross sections presented earlier 20 . As seen, the overall 
behavior of -RJ^f is well produced with fairly good fit qualities (see Table 5). In the case of the 
reduced cross section ov*^ only the bins M x G [0.28,2] and [2,4] GeV (the largest (3) are affected, 
resulting in a slightly better x 2 /dof. In Fig. 16 we show the corresponding total cross section alone: 
the match with the data is strikingly good indicating that the source of the uncertainties in is 
the incomplete description for the diffractive cross section. 



20 More R$j$ data with Q 2 e [2.7, 55] GeV 2 can be found in [73]. The fits to this data turned out to be slightly 
worse. 
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0.28 <M <2 



2<M <4 



4<M <: 



8<M <15 



15 <M <25 



25 <M <35 




100 200 100 200 100 200 100 200 100 200 

Fig. 15: The ratio of the diffractive versus the inclusive cross sections according to Eq. 4.16 as a function of 
y/s~/»p for different values of Q 2 and bins of diffractive mass Alx- The smallest values of f3 are located in the 
top right corner: from there, following the panels to the left or down increases /?. The ratios are based on the 
parameters for GT in Tables 1 (a to t) and 3 (crdiff , Eq. (4.10)). The data are from [41]. 



M x [GeV] 


0.28-2 


2-4 


4-8 


8-15 


15-25 


25-35 


E 


Treatment dof 


54 


57 


46 


37 


23 


10 


227 


if (3) ; N?(r): (4.5) 


2.05 


1.22 


2.27 


0.77 


1.16 


0.29 


1.51 


Ff<® ; N?(r): (4.10) 


2.05 


1.22 


2.30 


0.82 


0.88 


0.24 


1.49 


a? (3) ; Ny(r): (4.5) 


1.78 


1.16 


2.30 


0.77 


1.16 


0.29 


1.44 


a? (3) ; JV?(r): (4.10) 


1.78 


1.16 


2.33 


0.82 


0.88 


0.23 


1.42 



Table 5: The fit results of R^l for each M x bin (GT, diffractive cross section vs. reduced cross section). The 
data are from [41]. 
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x, . = 0.02 



• ZEUS 1999-2000 
GT 



I I i I i I i I i l^jl 

50 100 150 200 250 

1/2 
TP 

Fig. 16: The j*p cross section as a function of ^/s 7 * p = \J Q 2 (l/xb] — 1) for different values of Q 2 . The 
data are from [ ], i.e. it is the total cross section part of i?tot f data shown in Fig. 15. As long as Xbj = 
Q 2 /(Q 2 + s 7*p) J; 0.02 the theoretical predictions match excellently with the experiments. The breakdown at 
large Xbj becomes evident as one approaches the lower left corner with fixed ^/s 7 * p « 45 — 65 GeV. 

5 Asymptotic versus pre- asymptotic fits, a comparison 

As we have discussed at length for the total cross section in Sec. 3, an asymptotic fit without the 
energy conservation included is not feasible. This is caused by a strong fit tension arising between 
too large an evolution speed at small Y and too shallow a dipole correlator shape at large Y. 
Once energy conservation is included, both features improve towards what is needed to match the 
data: evolution speed is uniformly lowered and correlators become steeper overall. The upshot is 
a fit whose main ingredients are determined perturbatively - both shape and evolution speed in 
the asymptotic region are predominantly 21 determined by the nonlinear structure of the evolution 
equation and its kernel. As reported, it works flawlessly up to Q 2 of 1200 GeV 2 , i.e. the largest Q 2 
values available at HERA for x < 0.02. 

The fits performed in [45] by contrast include running coupling effects but no energy conservation 
correction. They obtain a fit that is almost as good as the GT fits described in this paper. To achieve 
this, they must move away from the asymptotic pseudo-scaling region and use pre-asymptotic 
features of an evolution trajectory. This allows them to simultaneously satisfy the speed and 
steepness requirements of the dipole correlators and improve the fit quality over what is possible in 
the asymptotic (pseudo-scaling) domain without the energy conservation correction included. The 
features of the correlators along the part of the evolution trajectory used in the fit differ drastically 
from those in the asymptotic case and a detailed comparison is in order. We will perform this 
comparison for both the total and diffractive cross sections (which [45] does not consider), to 
illustrate once more that, with present theoretical limitations, there is no hope to use the more 
differential diffractive cross section to differentiate between theoretical approaches. 

In doing so we have to deal with secondary differences of the fit procedures in our case and in [45] . 
Where we have argued for the use of current quark masses and have included a charm quark 

21 Aside from regulator effects on the coupling which are visible due to the R a values inherent to the kinematic 
properties of the HERA experiments. 
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contribution successfully using (optionally) an m/ -independent remapping of x i— > x c g — + q2 s ^ , 
[45] uses three quarks with TO/ = 140 MeV and an m/-dependent remapping of i h> x™" w = 

x ^ Q2 m/ , which, if used with our solutions degrades the x 2 considerably. The fit and parameters 
of [45] were done for x < 0.01 instead of the x < 0.02 we have used. For comparison we have 
to restrict ourselves to the same range. This reduces the Q 2 range of available data from Q 2 < 
1200 GeV 2 to Q 2 < 650 GeV 2 . Fortunately the separation scheme used to define the running 
coupling contribution is the same in both treatments. 

For the pre-asymptotic fit scenario, inclusion of the charm quark is virtually impossible without 
reworking most of the ingredients. To simplify the comparison we have also left out the charm 
quark in the asymptotic fits we compare to explicitly - for its inclusion see Sec. 3.3. 

We have attempted to expose the effects of the secondary differences in the fit procedures by 
comparing a set of fit scenarios that permute some of these ingredients. The fits below are tagged 
as BK^ hys , BKf 40 , BKf^ cft , [BKf^ e «] rc , and [BK^" a ] IC and defined as follows 

• BKp hys : asymptotic (pseudo-scaling), energy conservation correction included, three quarks, 
mass pattern to/ = {3, 5, 105} MeV. An x remapping induces no discernible differences. 



BK^ 40 : asymptotic (pseudo-scaling), energy conservation correction included, three quarks, 
mass pattern to/ = {140, 140, 140} MeV. 



BK^g cti : asymptotic (pseudo-scaling) , energy conservation correction included, three quarks, 
mass pattern to/ = {140,140,140} MeV, remapping of x according to x i-> x, 



GB-W _ 

eff ~~ 



Q 2 +imj 

X 7=T5 ■ 



• [BK^ 4 'q°"] 1c : same as previous, but with the energy conservation correction omitted. This 
attempts a fit on the asymptotic line the fit of [45] eventually merges onto. 

• [BK^4Q Ci '] rc , the fit developed in [45]: pre-asymptotic, NLO corrections restricted to run- 
ning coupling corrections (energy conservation correction excluded), three quarks, mass pat- 
tern to/ = {140,140,140} MeV, remapping of x according to i 4 X ^~ W = x — + Q^ i ■ 
Evolution starts from a GB-W like initial state Sy^Lj,^ 01 -. (r) = exp [— (rQ s . /2) 2 ] with 

Q 2 o = 0.241 GeV 2 . The scale choice for the running coupling in the fit is parametrized 
differently from our treatment: In the argument of a running coupling we use an r dependent 
scale in the form ^p- = with C 2 = 4e~ 27E ~3 and vary A. They use ^^p- = 4^, with 
A set to .241 GeV a priori, and obtain C 2 = 5.3 from the fit, see Fig. 17, right panel. 

In all treatments including [BK^|g off ] rc an APT regulator is used while [45] regulate the coupling 
by freezing it at a™ ax = .7, see Fig. 17, right panel. Note that the choice of initial condition in 
[BKj^o'^lrc is quite restrictive and no effort is made to vary the shape of the correlator other than 
allowing for some offset rapidity before matching evolution results to data. A systematic study of 
the impact of varying the shape, even a theoretical exploration of which shape features might be 
responsible for what kind of physics property of the cross section is still outstanding. Here we only 
attempt to contrast asymptotic fits with one example of a preasymptotic one. 
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Xbj < 0.01 




ni/A 
alV phys 


DT/ A 
BK -140 


BK 140 


[£5JV 140 J rc 


[BK 


Q 2 < 45 
GeV 2 


A [MeV] 
X 2 /200 
[GeV" 2 ] 


87.3 
0.82 
55.68 


59.0 
0.95 
72.3 


52.7 
1.01 
77.3 


31.2 
2.23 
79.3 


241 

0.97 

81.6 




A [MeV] 


93.7 


68.9 


63.1 


52.4 


241 


Q 2 < 650 
GeV 2 


X 2 /230 
[GeV^ 2 ] 


0.86 
55.91 


1.02 

69.9 


1.09 
74.2 


3.42 
71.4 


1.01 
81.7 



Table 6: A comparison of the fits to the total cross section. Parameters are based on the full range < -01, 
Q 2 < 650 MeV. The superscripts "A" and "P" refer to the asymptotic and pre-asymptotic fits, respectively. 
The quark masses are m u ,d,s = {3, 5, 105} MeV (physical) or m u ,d,s = {140, 140, 140} MeV (ad hoc). 



We intentionally only show BK based fits to allow for a direct comparison, despite the fact that 
GT fits have better x 2 and note that the modifications from BKp hys through BK^ 40 to [BK 
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incrementally reduce fit quality. Fit quality recovers only for [BK^|g off ] rc and relies on the freedom 
gained once one allows for correlator shapes away from the pseudo-scaling behavior. See Table 6. 
The first qualitative difference of the fit scenarios is captured in a plot of evolution speeds (left 



140 



and 



panel in Fig. 17) which indicates the parts of trajectories used in the fit of [45], [BK 
our favored asymptotic fits. The R S A ranges of the successful pre-asymptotic and asymptotic 
fits are strongly shifted against each other as required by the constraints on evolution speeds. 
Note that the fit interval matched onto data on the pre-asymptotic fit trajectory ends before the 
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Fig. 17: Left: Comparison of evolution speeds of pre-asymptotic fits without energy conservation correction [45] 
with the asymptotic fits including the energy conservation correction. The marked areas indicate the ranges of 
X(R S (Y)A) where the different fit strategies match to data. The curve starting from i? s ,oA. = 0.087 corresponds 
to the pre-asymptotic fit. Middle: The situation after the ratio C/A is known from the data fit. Due to the 
unique energy dependence of each scenario, the fit range 5 x 10~ r < Xbj < 0.01 is emphasized by thicker line 
width for each curve. Right a s (/j) for each scenario after the ratio C/A is determined by the data fit. C/A ~ 5.2, 
7.1 and 19.1 for the scenarios BK 



phys 



BKfVn and [BK^« 



c respectively. This gives an impression of the size of 
the slowdown effect originating from the energy conservation correction. The black dashed line shows coupling 
used in [ ], regulated in the IR by freezing it at maximum value of .7. 



asymptotic shape is reached (and the correlator shape would have become too shallow): the fit 
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strategy applied in the pre-asymptotic fit [BK^|Q 0ft ] rc [4- r >] relies on features of the initial condition 
throughout. Strikingly, the evolution speed increases over the whole fit interval of the pre-asymptotic 
fit, while it monotonically decreases in the asymptotic case. In addition, the dipole cross section 
of [BK^4Q Ctt ]rc interpolates between a steep, almost GB-W shape at Zbj = -01, and a final shape at 
the small Xbj end of the data range that remains still steeper than the (too shallow) pseudo-scaling 
shape obtained from running coupling BK evolution without the energy conservation correction . 
(See Fig. 18, left panel.) Note that detailed short distance information in this fit is predominantly 




r (scaled) r | GeV ' ] P 



Fig. 18: Correlator shapes obtained from the fits to the total cross section (Table 6). Left: A comparison of the 
shapes of iVy?(r). The curves are shifted to match at Ny q (r) = 0.5 for each value of Y. Middle: Ny? (r)/r 2 : 
the behavior at small r. Right: Contributions to the diffractive cross section for xp = 10 using a factorized 
b-profile. 



imprinted through the choice of initial condition. Evolution effects, which propagate into the UV 
just like in the linear BFKL case, have not yet reached the short distance tail of the correlators. 
This is different in our asymptotic fits, where the shape is determined entirely by the structure of 
r.h.s. of the evolution equation, which in its entirely is based on a (highly resummed) perturbative 
calculation. 22 To highlight the actual differences at short distances we show part of the short 

distance asymptotics of the fits by plotting ^ 2 in Fig. 18, right panel. 

We conclude our comparison with a look at the diffractive cross sections, see Table 7. The pattern 

ir P <0.01 B d a s /3<0.5 /3 > 0.5 Total Q 2 < 10 10<Q 2 <45 Q 2 >45 

X 2 /dof X 2 /153 x 2 / 21 5 X 2 /368 x 2 /97 X 2 /174 x 2 /97 



BKf 40 , fact. 6.30 0.40 1.39 1.37 1.38 1.53 1.35 1.4 



[BK[ 40 ] rc fact. 7.47 0.34 1.29 1.22 1.25 1.50 1.17 1.15 

Table 7: The fits to the diffractive data with the parameters shown in Table 6. Tag "fact." refers to the 
globally factorized case (4.5), however, with the large- N c replacement Cj — > N c /2. 



of fit quality shown in the table repeats what we have already seen for the total cross section: 

22 Operationally this is imprinted on the solution by having the solution evolve into the asymptotic region over a 
long Y interval before the fit range is reached and the correlators have reached pseudo-scaling shapes. Practically 
this will only affect the UV within the BFKL diffusion radius, but this easily covers the range relevant to HERA fits. 

23 For a direct comparison in momentum space, consult App. D.2, in particular Fig. 24, both the middle and right 
panels. 
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• The asymptotic fits shown in Sec. 3 and 4 show a good fit quality, also when restricted to 
£bj < 0.01 and Q 2 < 650 - i.e. the good fit quality shown earlier is not driven by the 
large x and Q 2 part of phase space. The choice m u ,d,s — 140 McV for the light flavors, 
although popular in the literature, is not optimal and the redefinition of the Bjorken variable 
x c g = Xbj(l + 4m 2 / Q 2 ) actually reduces the fit quality even further. Therefore, in Figs. 18 the 
curves corresponding to the modified BK evolution are based on BKp hys in which x e fi = Xbj- 
The fit to the diffractive data is fairly good. 

• Asymptotic fits without the energy conservation correction ([BK] rc of Sees. 3 and 4 or [BK^ 40 ] rc 
introduced here) result in a poor fit even within the range x^ < 0.01 and Q 2 < 45 GeV 2 . The 
fit within the wider range Q 2 < 650 GeV 2 is a failure. The modification m f = 140 — > 5 MeV 
improves the fit slightly, giving x 2 /200 = 1.65 with A = 52.7 MeV and a Q = 56.1 GeV -2 . 
The diffractive data are not considered. 

• For the pre-asymptotic procedure [BK^4Q° ft ] rc of [ ], a good fit is obtained for Xbj < 0.01 
and Q 2 < 45 GeV 2 . The fit quality remains good when the results are extrapolated to 
Q 2 = 650 GeV 2 . The fit to the diffractive data turns out to be excellent. The ratio a^/ap^p 
is not investigated. It is clear, however, that if both the diffractive and the inclusive cross 
section data are well resolved, then the same holds also for their ratio. 

Close inspection of Table 7 reveals that [BK^4g ot '] rc shows the best ^-vahie for diffractive data 
presented in this paper. The origin for this behavior can be discerned from the right panel of 
Fig. 18, which presents the diffractive structure functions obtained with the asymptotic and pre- 
asymptotic approaches. The overall shape of each component is the same in both approaches but 
the magnitudes are systematically smaller for pre-asymptotic case. This is especially pronounced 
at small (3 and explains the somewhat better match with the data. Unfortunately, it is the small 
j3 region where our present theoretical input, in particular the expression used to approximate the 
qqg contribution is deficient. Thus one cannot make any drastic conclusions based on this region: 
further corrections may change the results fundamentally. Our present theoretical setup is not firm 
enough to make use of the more detail information inherent in diffractive data to further constrain 
our analysis. 



6 Conclusions 

The main message from this study is that JIMWLK evolution at NLO allows an asymptotic fit to 
all HERA data below Xbj < -02, to both total and diffractive cross sections. This becomes possible 
only after all NLO corrections - including the energy conservation correction - are included. The 
energy conservation correction is the decisive ingredient in this argument: an asymptotic fit without 
it barely works for Q 2 < 45 GeV 2 and fails entirely beyond. Contrary to that, with the energy 
conservation correction included, the fit obtained for Q 2 < 45 GeV 2 simply extrapolates to a 
successful fit all the way up to Q 2 = 1200 GeV 2 , i.e. the low Q 2 range determines the result 
fully, over a range much larger that expected from naive BFKL based momentum space diffusion 
arguments. While our study does not exclude that a pre-asymptotic component is compatible with 
data also if one includes the energy conservation correction, such a feature is not required with 
present accuracy. 
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This observation is quite striking in that such an asymptotic fit is highly constrained and largely 
determined by perturbation theory: the shape of the asymptotic, pseudo-scaling correlators them- 
selves is fully determined by the structure of the r.h.s. of the evolution equation, which is the result 
of a purely perturbative calculation. 

The pre-asymptotic fits of [45,46], which omit the energy conservation corrections, on the other 
hand, rely strongly on the fact that the evolution trajectory matched onto data has not yet reached 
the asymptotic domain. Both evolution speed and correlator shapes can only be matched onto data 
as long as many of the features of the initial conditions are not yet erased. This remains true even 
though the precise details of these initial conditions are not too strongly constrained as soon as one 
allows evolution trajectories that start with shapes not related to the GB-W model by evolution. 

This puts an additional emphasis on the role of NLO corrections: without at least the inclusion of 
running coupling corrections a fit to HERA data from evolution equations is virtually impossible - 
it becomes very hard to find an initial condition that would allow us to match the full HERA range, 
mostly due to too fast evolution. NLO corrections bring both qualitatively new features (such as 
scale breaking via the running of the coupling) as well as quantitative modifications (slowdown 
via the energy conservation correction) that make a data fit successively easier, until, with both 
included, even the UV details are naturally imprinted by the evolution equation itself. 

All other results are secondary to these observations, but they round out the picture and point us 
to what main theoretical improvements are needed to step beyond what can be done presently. 

The total cross section fits work best if one assumes massless quarks - a fit with physical quark 
masses starts to deviate from at Q 2 values below mass of the heaviest quark included in the fit. 
Phenomenologically only u,d,s, and c quarks need to be considered and only the charm quark 
with a mass of the order of Q s induces a strong modification. Since all these modifications are 
concentrated in the infrared, one should either drop this region from consideration altogether and 
accept the loss of fit quality in this region, or adopt non-perturbative arguments to improve the fit. 
One such strategy improves fit quality in the infrared by distorting how one maps the evolution 
trajectories onto phase space by replacing x^ by some <5 2 -dependent x e g. Even with rough models 
for x e g, plausible charm fractions are achieved. Any refinement that aims at quantitative rather than 
qualitative improvements, however, must incorporate nonpcrturbative information from outside the 
scope of JIMWLK evolution. 

The asymptotic fits to diffractive data use the parameters already determined via the total cross 
section. The formalism clearly requires contributions from NLO impact factors to fill in phase 
space at f3 — > 0. Without them a fit to data is not possible. While this was known empirically 
already from fits within the GB-W model, from an evolution perspective this becomes a consistency 
requirement. As soon as more than one Fock- component comes into play the details of how one 
models the impact parameter dependence of their scattering begins to affect fit quality: the choice of 
model for this non-perturbative aspect starts to affect the relative weight of these Fock-components 
in the associated cross section. One of the few direct experimental constraints in this respect comes 
from the ratio of total to diffractive cross sections, and can be accommodated easily with Gaussian 
profiles no matter if one uses a globally factorized form or enforces exact Casimir scaling on the 
other extreme. 

At present, two issues hamper precision fits that involve differential cross sections with exclusive 
final states: the lack of complete NLO impact factors, and a consistent treatment of the impact- 
parameter dependence based on non-perturbative input. Only the first of these issues is sure to be 
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resolved in the near future, but it requires considerable effort: The full NLO corrections to JIMWLK 
evolution and the NLO impact factors have to be calculated. Beyond that a treatment of resummed 
collinear corrections needs to be formulated that eliminates any double counting issues that affect 
our current inclusion of the energy conservation corrections. All of these can be addressed within 
perturbation theory. Any progress on the impact parameter dependence requires is a different 
matter altogether. 
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A Kinematics and common approximations 

For completeness we include a brief description of the kinematical variables involved in the inclusive 
and diffractive deep inelastic scattering processes. 

Only two independent variables are required for the process 7* (q) p(P) —> anything, the photon 
virtuality Q 2 and the Bjorken variable x: 

Q 2 := -q 2 := -{k - k') 2 (A.la) 

where, at leading twist, Xbj carries the interpretation of the momentum fraction carried by the quark 
inside the target that is struck by the virtual photon. The reduction to two kinematic variables 
holds, wherever 

P 2 /Q 2 = m 2 p /Q 2 < 1 and l/z bj > 1 (A.2) 

to guarantee that s 7 * p , the total energy squared of the 7* (5) p(P)-subprocess (frequently denoted 
W 2 ) can be expressed in terms of Q 2 and x^j only: 

s 7 * P = (P + q) 2 = P 2 + 2P-q-Q 2 ^Q 2 ( — - 1 ) . (A.3) 

The conditions (A.2) are satisfied for the majority of the experimental data. 

In addition to the variables in Eqs. (A.l), two additional kinematic variables are required to describe 
the kinematics of the diffractive process 7* (5) p(P) X{Mx) p{P')'- 

P P-q Q 2 + STp -m 2 p ~ Q 2 + s Tp ' lA - 4aJ 

R= Zl = 91 Q 2 = Xbi (A4b) 

P 2(P-P')-q Q 2 +M 2 x -t Q 2 +M 2 X x P ' { ' ' 
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Fig. 19: x 2 /dof as a function of A(Y" ff). Each evaluated x 2 /dof 
corresponds to a different choice of N(Y — 0,r) with the op- 
timal A and <tq. The line at x 2 /dof = 1 shows the extraction 
point of the error bars for A (see Figs. 7). In the asymptotic 
approach the offsets Y ff are actually negative meaning that one 
discards a certain number of configurations from the beginning 
of the evolution (typically Y = is in the pseudo-scaling stage). 
The horizontal crossing line (blue) indicates the normalization 
dependence via 0"o/o"o,opt wnere o"o,o P t = 55.33 GeV~ 2 is the nor- 
malization of the best overall x 2 /dof. This log-linear dependence 
of the normalization is common to all descriptions. Note the red 
dot-dashed curve in the top left corner that corresponds to the 
data fit without the energy conservation correction. 



where xp is the fraction of the proton four momentum carried by the colorless diffractive exchange 
called the Pomeron. The variable (3 is an analogue of the Xbj for the diffractive system: it is 
the momentum fraction of the Pomeron carried by the interacting parton inside the Pomeron. 
The variable Mx is the invariant mass of the diffractive hadronic final state denoted by X and 
t = (P — P') 2 < is the squared four momentum transfer. Experimentally \t\ <§; Q 2 , M x , thus t 
is set to zero in the above equations in addition to the proton mass m p . 



B Parameter optimization 

The total error for the inclusive data is obtained by adding the systematic error quadratically to 
the statistical error as follows, 

etot = ±J(stat)2+max[(sys ± )2] ; X 2 /dof = ^ ~ /dof , (B.l) 

i \ e tot) 1 

which results in symmetric error bars. For the inclusive data this is fine since only a small number 
of data points have asymmetric systematic errors. However, in the case of the diffractive data the 
asymmetric error bars are kept due to the large differences between |sys + | and |sys_|: a bias of 
this size in the error propagation cannot be ignored. A common % 2 -test for the goodness of fit is 
applied, shown on the right in (B.l). The parameter optimizing procedure is illustrated in Fig. 19 
where x 2 /dof is plotted as a function of the (correlated) parameters A(l^ff). 

In some figures and tables we show the error bars for the BK based descriptions. The origin of 
these is the data fits and they are based on the arbitrary choice made by us. Even though there 
are three parameters to be fitted, it turned out that the procedure of seeking the minima of x 2 is 
effectively one parameter fit: two of the parameters are correlated and the third one is just a trivial 
normalization. We illustrate this in Fig. 19 where x 2 /dof is plotted as a function of the correlated 
parameters A(Y s). The horizontal dotted line indicates the point from where the error bars are 
extracted: it is chosen to be at x 2 /dof = 1. 
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C Numerical implementation of the energy conservation cor- 
rection 



The numerical treatment of the BK- and Gaussian truncations of JIMWLK evolution is identical, 
an appropriate invertible map relates the two truncations for qq-dipoles [51], despite any differences 
encountered in the application to more general correlators. This statement in fact extends to the 
inclusion of the energy conservation corrections as shown in Eqs. (2.32) and (2.33). Since the actual 
numerical implementation used is done in BK-form, we phrase the discussion here in BK form only. 
With the exception of the energy conservation correction, the discretization and the numerical 
evaluation of the evolution equation follow [51]. Since the computational cost of a single right-hand 
side evaluation scales like N paTcnt x A^aughter x Wangle (where Ni denote the number of grid points 
in polar coordinates), optimization of the number of discretization points is highly recommended 
when solving the modified BK equation 24 . 

The F-derivative appearing on the r.h.s. of Eq. (2.32) can not be evaluated directly, to cope with 
it we have implemented an iterative procedure. We first write the evolution equation in finite 
difference form 



ck _ q. 
SY 



f[Si. 



/fe 1 ] - f[Si. 
SY 



F[St?] ; f[S] — r.h.s. of BK , (C.l) 



where the second term on the r.h.s. represents the energy conservation correction. 

In (C.l), S is a function of dipole size and i labels discrete rapidity steps, k is used to label iteration 
steps at fixed rapidity used to deform the solution of the equation without energy conservation 
correction, i.e. S® +1 , defined via 

%^ = f[Si] (C.2) 

into a solution of the full equation (C.l) at Si+i := Sk^°°, where the limit k — > oo assumes 
convergence of the procedure. Leaving convergence issues aside for the moment, the iteration 
with (C.l) proceeds as follows: 

1) evaluate f[S%£] 2) calculate F[S^] 3) evolve Sj +1 = Sf + SYF[S^} 

4) return to 1) with new configuration from step 3). (C.3) 

The iterative solution is accepted as the solution Si+i for the next time step at some finite k, once 
left and right hand sides of the discretized equation (C.l) agree to some desired accuracy e 

qk _ qO 

i+1 sy * -F[S? +1 ]<e, (C.4) 
or equivalents F[S* +1 ] - F[S^} < e. 

24 High computational cost of the R.H.S evaluation in the iterative procedure together with the requirement a small 
rapidity step 8Y makes the overall numerical cost of a single evolution trajectory far too large. 
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Note that the definition of evolution speed A in Eq. (2.39) can be used to directly translate the 
accuracy criterion of Eq. (C.4) into the error implied for A: 



A f := 2 



J ^[l.h.s(C.l)-r.h.s(C.l)] = 2 J 



dr 
r 



ok 



s? 



SY 



(C.5) 



The main issue with this procedure is that it is not stable for arbitrary shapes of the dipole function 
S: convergence towards (C.4) occurs only for a small number of iteration steps before the iterations 
start to diverge in the sense of an asymptotic series. Convergence can be improved by reducing 
SY (to, in general, impractical values) with the pseudo-scaling region showing the best convergence 
properties. 

A numerical exploration reveals that it is the nonlinear region r 2 > i? 2 (Y"), where Ny(r) = 1— SV(r) 
approaches one, that is least stable. The source of the instability is the nonlinearity, but the 
convergence properties of the asymptotic series and the obtainable accuracy can be improved by 
reducing the step size SY. 

The manner in which SY affects convergence of the iteration procedure is slightly peculiar since 
the energy conservation term [in the update step (Eq. (C.3), step 3))] carries no explicit overall 
power of SY. However, its initial size (in the first iteration step at k = 1) is determined by the 
difference between Si and S® +1 as induced by (C.2). This difference is proportional to SY and thus 
step size imprints itself on the whole subsequent iteration procedure: convergence can be improved 
by reducing SY. 

The main criterion for convergence therefore is the difference between Si and Sf +1 , and this not 
only depends on SY itself, but also on the shape of Sf. a solution near the pseudo-scaling regime 
generically leads to a smaller energy conservation correction and better stability than a solution 
that is far from pseudo-scaling, such as the Gaussian shape used by GB-W. 

Since stability is most precarious at r 2 > i? 2 (Y"), there is an interplay between IR regulators applied 
to tame the Landau pole and the convergence of the iteration procedure: Irrespective of the shape 
of the solutions, stability is improved whenever SYK xzy R c xzy is small. R x s zy plays the role of an 
effective coupling and we find that 



• the linear region r < R S (Y) is stable even with the relatively large fixed coupling a s = 0.4. 
A fast convergence occurs in all cases. 

• the non-linear region r > R S (Y) is unstable even with the relatively small fixed coupling 
a s — 0.2. The case a s = 0.4 is already challenging and requires an impractically small step 
size SY. The convergence is generally slow. 

For the realistic case with the full running coupling kernel (whose size is essentially determined 
by the size of the parent dipole r) this implies that the problems caused by the presence of the 
nonlinearities are exacerbated by the running of the coupling. Thus, the iteration is quick and 
stable at r < R S (Y) thanks to both the smallness of the kernel and absence of nonlinearities. In 
contrast to that, the region r > R S (Y) is difficult (especially at low rapidities) since the kernel is 
large and the equation is dominated by non-linear effects. 

In practical terms, the convergence properties of the iteration procedure preclude the use of (C.3) 
away from the asymptotic line. We have applied (C.3) in the pseudo-scaling region, iteratively 
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reducing SY to push the step of minimal error fc m i n to larger k and minimize the error in a brute 
force approach to set a baseline. As is typical with iteration procedures, convergence properties can 
be strongly affected by a modification of the iteration procedure. As a compromise between speed 
and accuracy near the pseudo-scaling region we have amended steps 1) — 2) of (C.3) by a set of 
re-weighting steps: 

o\ ok qO i sy rro^ — 1] /~ik O — 1 Qk 

4) evaluate /[St+i] and calculate F[Si +1 ] 

5) S& 1 = S° + 5YF{S* +1 ] ; + 1 = 2S* +l 

6) calculate Wf +1 = - A s (G* + 1 - / (A, - A.) 
where A, = £(S#-i- ; A s = ^G^ 1 - 

r r 

7) return to 1) with weighted W^ +1 of 6) taking the place of (C.6) 

The procedure terminates when||A s /A g | — 1| < e. Then further iterations of do not lead to 

any improvement 25 and one proceeds with Sf +1 in step 4). Initially (typically) |A S | ^> |A S | and so 
Gf +1 gets weighted less than Sf +1 in W^ +1 . 

The upshot is high precision at small r after few iterations already at SY = 0.0025, at the price of 
comparatively limited precision in the large r part and with very few iteration steps. The procedure 
always terminates but can, without additional modifications, not exceed the precision indicated in 
Figs. 20 and 21. 

We observe that the unmodified method (C.3) at k = 1, i.e. the crudest approximation, greatly 
underestimates A whereas the first iteration k — 2 overestimates it, as seen in the left panel of 
Fig. 20. 

From k — 3, the brute force solution to (C.3) is more accurate at large r than the approximative 
solution (C.6), however, this solution is increasingly difficult to obtain: in order to keep the iteration 
stable one has to reduce the change in shape of S with Y by working near the pseudo-scaling region, 26 
and in addition needs to employ a very small rapidity step SY < 10 -3 . To stabilize (C.3) also for 
k = 4, 5 requires SY < 10~ 4 ,10~ 5 respectively which increases the numerical effort prohibitively. 
This is illustrated in the right panel of Fig. 20 which compares the three step iteration of (C.3) with 
the result of the re-weighted result from (C.6) obtained with almost an order of magnitude smaller 
CPU time. Fig. 21 shows the size of the energy conservation correction in the different cases. For 
both Figs. 20 and 21, the configurations are selected from the Y- or i? s A-range typically used in 
the data fits 27 . 

The induced error of evolution speed according to (C.5) is shown in Table 8 for three samples 
corresponding to R S A — {0.3,0.2,0.1}. The region r < R s is fairly accurate in all cases and 
especially the case k — 3 is close to the accurate solution even at r > R s . The peak values 
(Figs. 20, middle) correspond to Sy, X y ~ 0.2. As seen, the iterations k — 4, 5 do not bring any 

25 The condition A s pa — A g yields W^ +1 ~ + G* +1 )/2 regardless of |A S>9 | <S 1 which obviously is a bad 

solution for any r. 

26 This requires some numerical optimization of its own. 

27 In the data fits A w 0.1 GeV and so R S A = {0.1,0.2,0.3} -¥ R s w {1,2,3} GeV -1 . 
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Fig. 20: Left: A: a comparison of iterations k = 1,2,3,5 of (C.3), the average of k = 1,2 and the proce- 
dure (C.6). Middle: Accuracy of the solution from (C.4) as a function of r for R S A. = {0.1, 0.2, 0.3} (corresponds 
R s ~ {1,2,3} GeV - ). Right: The solution of the case k = 5 is, for all practical purposes, accurate but the 
requirement of computing time is huge. L.H.S. and R.H.S. refer to the discretized evolution equation (C.l) 

SY A e |A e (r < R s )\/X e 

reweighted 0.0025 {0.025,0.018,0.011} {0.256,0.250,0.245} 

brute force k = 3 10~ 3 {0.010, 0.007, 0.003} {0.106, 0.092, 0.079} 

brute force k = 5 10~ 5 {5.824, 3.541, 1.561} x 10~ 3 {6.437, 0.448, 4.590} x 10~ 3 

Table 8: The error estimates for A based on (C.5) for different iteration procedures (first column). The 
numerical cost is measured by SY. Most of the error comes from \r\ > Rs(Y) (last column). The values 
correspond to correlation lengths R S A — {0.3,0.2,0.1}. 



essential improvement and for k > 5 the iteration procedure turned out to be highly unstable at 
r > R s for any relevant SY. 

It should denoted that a convergent solution in the offset region 29 is easier to obtain with some 
other, say stronger, regulator but, however, inside the actual fit range this kind of modification does 
not bring any improvement. 



D Tools to efficiently address diffractive cross sections 

A lot of simplifications used in earlier treatments with simple height one box profiles can be carried 
over to more general profile models if the r and r' dependence in the b-integrated product of dipolc 
amplitudes factorizes according to 

J d 2 b N^ Y (r,b)N n , Y (r',b) = o% [N^ Y (r)}* N^ Y (r') , (D.l) 

or if such factorizcd behavior is a good approximation to the full result c.f. (4.10). 

28 The iteration seems to be approaching (alternatingly) the fixed point but the requirement of a very small SY 
makes the brute force method impractical. 

29 The region of large running coupling that is cut off by the parameter Y g . 

30 In (D.l) the constants cr^' are adjusted such that N^ Y (r) 1 > °°> 1 for convenience. 
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Fig. 21: The energy conservation correction (f[Sf + i] — 
f[Sf])/5Y for the cases k = 1,2,3 of (C.3), an average of 
k = 1,2 and the procedure (C.6). The latter approximately 
coincides with the average of k = 1,2 for all r and thus the 
easiest way to get a fairly good solution is taking the average 
of k = 1,2. A sharp slowly converging peak at large r is a 
consequence of the non-linearity since it appears even with a 
relatively small fixed coupling a s = 0.2. The running coupling, 
being large at these scales, gives rise to the additional stability 
problems. 



Wherever (D.l) holds, a momentum space variant of the cross sections of (4.1) offers an efficient 
route to perform the data fits. It allows to pre-calculate most of the numerical integrals in the 
final expressions for the diffractive structure function and speeds up the parameter seeking process 
tremendously. The tools to do so are collected in the remainder of this section. 
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D.l Diffractive cross sections in momentum space 

The momentum space expressions for the three terms in Eq. (4.1) read 



d/3 



x d 2 b 



(27r) 2 a e 



E< 



Q 2 (z 2 + z 2 )N c 



32 



do 



2 Mf(q,b) [1-2(3- 



2mj 

zzQ 2 



(l-2(3)^-(q 2 + 2m 2 ) 



( £ T+9 2 ) 2 -4(/3 



m 2 ) q 2 



(27T)V 

Q 2 P 



E 



dz 



(3m 2 N c 
8irzz 



d 2 b 



zzQ 2 



zzQ 2 



mV) q 2 , 



(D.2a) 



dp 



da 



d 2 b 



qqg,T 



(2n) 2 a c 
Q 2 P 

dq 2 



E< 



/34Q 2 (zz) 2 jV c 

&7TZZ 



zzQ 2 



q 



2 N™{qMl- 



(^ + g 2 ) 2 -4(/3 



TOf ) q 2 , 



(D.2b) 



d/3 
x ld 2 b 



(2n) 2 a em a s pC f N c 



LLQ 2 
d 9 2 



Af^(q,b)[z 2 + z 2 + ^- 



{q 2 -fc 2 (l-2z)) 2 + 2zzk 4 
k 2 ^(q 2 + k 2 ) 2 -4(l-z)q 2 k 2 



(D.2c) 



Where we have used the shorthand expressions z = 1 - z and zf = |(1 ± ^l-Amj/M 2 -). 

In the massless limit, these expressions match up with their counterparts in [31,69] once the 
corresponding fo-profiles 32 have been inserted and the impact parameter integral has been carried 
out. 

The cumbersome expressions in brackets arise from integrations over the orientation of q of rather 
simple Fourier expressions for the McDonald Ki appearing in the coordinate space variants. The 
solutions of the evolution equations enter via JVy{q, b), which are determined from the coordinate 



L Use a s J r (x, q) 



f^Nqq^xxti/x) (q) to translate the qq expressions into those of [31,69]. The qqg-teim in the 



original literature suffers from an additional incorrect rescaling by a factor of (Cf /N c ) 2 and otherwise substitute Nqq 
for M gg [53,71]. 

32 Where ever these papers make use of the da / dt\t—a in their expressions one must use factorized Gaussian profiles 
for this procedure, the corresponding expressions are valid only in this case. 
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Fig. 22: A comparison between the momentum 
and coordinate space formulae of the diffractive fi- 
nal states. Here, the same dipole input (a quark 
dipole) is used in all contributions and the param- 
eters Bd and a s are set to one. In the case of 
the coordinate space equations, three light flavors 
with equal masses m Ut d.s = 5, 140 MeV are con- 
sidered whereas for the momentum space equa- 
tions mf — (coincides perfectly with the case 
m u ,d,s = 5 MeV). The curves are calculated with 
fixed Q 2 = 15 GeV 2 and x v = 0.001. 



space dipole amplitudes Ny{r, b). These are not related by a direct Fourier transform, see Sec. D.2 
for definitions and properties. 

The benefit of using the momentum space forms of the equations shows up in the numerical imple- 
mentation: the required rapidity range of Af can be pre-calculated whereas the coordinate space 
expressions have highly oscillatory integrands in which no part of the nested integrals can be pre- 
calculated. 

We have cross-checked our results using both coordinate and momentum space variants numerically. 
A direct comparison of the quark contributions at different masses gives an idea of how they affect 
the cross sections: non-zero quark masses reduce x^F^^j whereas xpF®^ is practically unaffected. 
In any case, it is expected that considering the non-zero quark masses would only lead to a small 
rescaling of the normalizations in Eq. (4.1) 



D.2 Integral transformations for dipole amplitudes 

The key ingredient is the non-standard "Fourier" -transform of the dipole amplitude (here adapted 
to the azimuthally symmetric forward case relevant in conjunction with (4.10)) 

N?(r) = J ^(l-e^)^( 9 ) (D.3) 

(1Z = qq, gg etc. labels the representation; see [62] for a broader exposition on the structure of the 
exponentials) together with its inverse 33 

N?{q) = (^2)) W(<7) where fa Y (q) := ± J ^V^iV£(r) . (D.4) 

Note that N as used here is dimensionless just as the dipole amplitudes. This differs from the 
analogous quantities in [(>:] or in the work of GB-W. We note that Af/q 2 is normalized to the 

33 This can be derived via an inverse Mellin transform [74,75] (see [76] for a step by step exposition). 
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saturation value of the dipolc amplitude: 

/ ^V?(9) = -Ny(M->oo) ; (D.5) 

correspondingly, Sy(r) = J" c ^e irq Ny(q). Eq. (D.5) provides a stringent check for our numerical 
tools, which also faithfully resolve the chain of transformations Ny(r) — »• 4>y(q) Ny(q) f° r the 
GB-W model, where all the steps can be determined analytically from Eq. (D.4). 

The chain of transformations Nyir) — > (fry (q) — > Afyio) is illustrated in Fig. 23, quark and gluon 
dipoles are compared in Fig. 24. 





r | GcV"' ] q [ GcV ] q 1 GcV ] 

Fig. 24: A comparison of the quark and gluon dipole inputs. From left to right: Ny(r) — > 4>y{q) ~^^y(q)- 
Green dot-dashed curves are for the pre-asymptotic fit scenario studied in Sec. 5 (starting from Y = ln(l/0.01) ss 
4.61). 
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E Consistency checks 



With the importance of NLO contributions firmly established, one should at least attempt to 
understand in which sense the results obtained from (2.32) or (2.33) are stable against modifications. 
One obvious modification is a simple numerical check for stability against higher order corrections 
such as higher order running coupling contributions. While this has already been discussed in 
the quite sophisticated framework of renormalon corrections in [54] a brief numerical check on the 
quantitative impact of such corrections on a data fit gives an alternative ballpark impression of 
their impact. 
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Fig. 25: The slowdown due to the running coupling is stable against higher order corrections. The vertical lines 
bracket the i? s -range of HERA. Left: One loop and two loop running couplings coincide after readjusting A. 
Middle: Evolution speeds differ only slightly in the HERA arrange. Right: The dipole correlators (magnified 
by r~ 2 ) agree after adjusting the fit parameters with only tiny differences remaining at short distances. 



The fit is repeated by using the two loop running coupling 

1 ft ln(ln( M 2 /A 2 )) 



c? (M) = £ 
Po 



ln( M VA2) Pi lnV/A 2 ) 



ft = 9, ft = 64 for N f = 3 , (E.l) 



which was also regulated by the APT method 34 . Once the optimal Ail^l and the corresponding 
rapidity offsets are known, the one and two loop couplings approximately coincide (see Fig. 25, 
left). Furthermore, the evolution speeds (Fig. 25, middle) and the dipole correlators (Fig. 25, right) 
coincide around Y = 10, which is a direct consequence of the fact that most data points are located 
around this rapidity. The fit quality is not essentially affected, see Table 9. 

The sensitivity to the infrared regulator is investigated. Whereas the APT regulator ((E.2), right) 
offers a smooth crossing over the Landau pole with the limit af PT (0) = 4tt//3o, for instance, adding 
of a constant inside the logarithm as follows, 

alt _ Am 1 apt/ \ _ / 1 1 \ , . 

a ° W ~ ^ln(^/A 2 + #) ; Us ^-Mln(/* 2 /A 2 )~M7A 2 -lJ ' { ' 

results in a steeply increasing (for not too large jf) coupling near the Landau pole. As seen in 
Fig. 26, the effect induced by this type of modification extends to the actual fit range. Despite the 

34 In practice, a 2L {^i) must be evaluated by using a spectral integral representation, see [63,64]. 
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GT ; a; bj < 0.02 


Q 2 < dof 


X 2 /dof 




[MeV] 


<x [GeV" 2 ] 


A 2 l/Ail 


One loop 


45 GeV 2 224 


0.81 


104.7 




55.33 




One loop 


1200 GeV 2 295 


0.80 


104.2 




55.24 




Two loop 


45 GeV 2 224 


0.81 


213.8 




55.41 


2.04 


Two loop 


1200 GeV 2 295 


0.83 


229.6 




54.59 


z.zu 


Table 9: Fit results for GT with two loop running coupling. The quark 


masses are m u 


, diS = 5 MeV 


BK ; x hi < 0.02 Q 2 < 


dof 


X 2 /dof 


A [MeV] cr [GeV" 


"I 


regulator: 


45 GeV 2 


224 


0.81 


93.7 


55.05 




APT 


1200 GeV 2 


295 


0.80 


97.7 


54.50 




regulator: 


45 GeV 2 


224 


0.96 


111.7 


56.57 




alt. 


1200 GeV 2 


295 


1.01 


102.9 


57.24 





Table 10: Fit results corresponding the regulators in Eq. (E.2). The case of the alternative regulator corresponds 
to the choice # = exp(l) with the limit af*(0) = 4tt//? . 



APT 

# = exp(l) 

# = exp(2) - 1 

# = exp(3) - 1 



a = 4lt / ( P ln( |i 7A~ + # ) ) 
R A = 2 exp( -5/6 -y ) 



it 
■ I 




Fig. 26: A comparison of the alternative regulators and APT 
regulator presented in Eq. (E.2). The effect on the evolution 
speed is not what one would naively expect it to be: on the 
one hand, a weaker regulator with a larger effective running cou- 
pling i?. e fF produces a faster evolution speed at the pre-asymptotic 
stage than the APT (as it should) but suddenly drops below it 
inside the fit range R S A ~ 0.1 — 0.4. On the other hand, a 
stronger regulator, and hence a smaller i? e ff, extends the pre- 
asymptotic stage remarkably which eventually leads to a faster 
evolution speed inside the fit range. The evolution with the 
APT regulator matches the best with the data. The length scale 
R a A = 2e _5//6 ,E indicates the location of the Landau pole. 



deviating evolution speeds, all cases shown can resolve the data with a good x 2 /dof < 1. However, 
as seen in Table 10, the evolution with the APT regulator yields the best fit to the data. 

The last check concerns the solution of the modified BK equation introduced in App. C. The fit 
with m Ut d, s = 5 MeV is repeated by using more accurate solutions for the modified BK equation, 
i.e. the cases k = 3, 5 of the procedure (C.3). The results of these fits are presented in Table 11. 
In both cases x 2 /dof is increased if compared with the reweighted case shown in the first row but, 



BK ; Xbj < 0.02 xV 2 95 A [MeV] rj [GeV -2 ] Table 11: Fit results corresponding to 

the different iterative solutions presented 
reweighted, Eq. (C.6) 0.80 97.7 55.50 in App c , n aN caseS| a wider data range 

k = 3, Eq. (C.3) 0.86 98.2 53.89 < i 20 GeV 2 is considered. 

k = 5, Eq. (C.3) 0.95 95.0 54.18 
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however, remains below one. The parameters A and <7n are altered as well but are still in the same 
ballpark as the ones shown in the first row. The experimental data seem to favor slower evolution 
speed obtained by the approximate solution. 



Reweighted 
k = 5 

Reweighted: rA -> r 
k = 5:rA->r 
N(Y,r) = 0.5 



" / / 




< / 






/ 




1 ©' / 




/ II / Fil 








/ */ / 


i 
/ 


— Reweighted 




— k = 5 


j 


— Reweighted: rA -> r 




--- k = 5: rA->r 


i 



— Reweighted 


rA 


— k = 5: rA 




— Reweighted 


rA -> r 


— k = 5: rA -> 





Fig. 27: Left and Middle: A comparison of Ny q (r) of the applied and accurate solution. For a fixed R 3 A, 
there is no visible difference between the shapes of the correlators. The same feature is approximately preserved 
after the physical scales are determined by the data fits. The correlators are extracted from R S A = {0.1,0.3}. 
The crossing blue dotted line is indicating the saturation condition Ny q (r) = 0.5. Right: The deviation between 
the corresponding evolution speeds remains after readjusting A. Thus, the reason for the deviation in x 2 /dof is 
a slightly different energy dependencies of the saturation scales rather than the actual shapes of the correlators. 



To summarize, the fits based on the modified BK/GT evolution are stable against a large variety 
of modifications. As seen in Fig. 26, the biggest uncertainty to evolution clearly emerges from the 
infrared regulator. 
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